🔢 mathematics

Subspace gradient descent method for linear tensor equations

This paper proposes two new subspace gradient descent methods for solving symmetric positive definite linear tensor equations using the Tucker format, mixed-precision strategies, and preconditioning, demonstrating their competitiveness against the state-of-the-art AMEn algorithm through experiments on 3D partial differential equations.

Original authors: Martina Iannacito, Lorenzo Piccinini, Valeria Simoncini

Published 2026-02-26

📖 4 min read🧠 Deep dive

CC BY 4.0

Original authors: Martina Iannacito, Lorenzo Piccinini, Valeria Simoncini

Original paper licensed under CC BY 4.0 (http://creativecommons.org/licenses/by/4.0/). ✨ This is an AI-generated explanation of the paper below. It is not written or endorsed by the authors. For technical accuracy, refer to the original paper. Read full disclaimer

The Big Picture: Solving a Giant 3D Puzzle

Imagine you are trying to solve a massive, multi-dimensional puzzle. In the world of math and science, these puzzles often represent real-world problems like predicting how heat spreads through a 3D object, simulating chemical reactions, or modeling financial markets.

Mathematically, these problems are written as Tensor Equations. If a standard math problem is a flat sheet of paper (2D), a tensor is a giant, multi-layered cake or a 3D Rubik's cube. The goal is to find the "perfect shape" (the solution) that fits perfectly into the puzzle.

The problem? These puzzles are so huge that they would crash a normal computer if you tried to store them all at once. It's like trying to carry a library of books in your backpack; you need a smarter way to pack them.

The Solution: The "Tucker" Suitcase

The authors of this paper propose a new way to solve these puzzles. They use a method called the Tucker format.

Think of the Tucker format as a collapsible suitcase. Instead of carrying every single book (every single number in the data), you only carry the essential "core" items and a few "sleeves" (matrices) that tell you how to expand them back out when needed. This keeps the memory usage low, allowing the computer to handle these giant 3D problems without running out of space.

The New Strategy: Gradient Descent with a Twist

To solve the puzzle, the authors use a strategy called Gradient Descent. Imagine you are standing on a foggy mountain and want to get to the bottom (the solution). You can't see the whole path, so you take small steps in the direction that feels like it's going down.

The Old Way: Usually, you just take a step down the steepest slope. This works, but it can be slow and zigzaggy, like a hiker taking a very winding path.
The New Way (Subspace Gradient Descent): The authors created two new "hiking guides":
1. The Steepest Descent Guide: A fast, simple guide that takes direct steps down the hill. It's quick to start but might get stuck in a loop if the hill is tricky.
2. The Conjugate Gradient Guide: A smarter guide that remembers where it has been. It doesn't just go down; it calculates the best angle to avoid backtracking. It's like a hiker who uses a map to cut across the valley instead of walking around the mountain.

The paper shows how to make these guides work specifically for those giant 3D "suitcases" (tensors) without breaking them open.

The Secret Weapons: Preconditioning and Mixed Precision

To make these guides even faster, the authors added two "power-ups":

Preconditioning (The GPS): Sometimes the mountain is so steep or the fog is so thick that even the best hiker gets lost. Preconditioning is like giving the hiker a GPS. It doesn't solve the puzzle for them, but it reshapes the mountain so the path to the bottom is straight and easy. The authors tested different "GPS" systems (like using Fast Fourier Transforms, which are like magic math shortcuts) and found that some made the solution happen in just a few steps instead of hundreds.
Mixed Precision (The Fast-Forward Button): Doing every single calculation with extreme precision is like measuring a mountain with a microscope—it's accurate but takes forever. The authors realized they could use a "ruler" (lower precision) for the rough steps and only switch to the "microscope" (high precision) for the final touches. This saves a massive amount of time and computer power.

The Results: Beating the Competition

The authors tested their new methods against the current "champion" solver, called AMEn.

The Test: They simulated heat flow and diffusion in 3D space (like heat spreading through a metal block).
The Outcome: Their new methods were faster and more efficient. In many cases, their "GPS" (preconditioning) allowed them to solve the problem in just 2 or 3 steps, whereas the old champion took many more.
The Trade-off: The "Steepest Descent" guide was sometimes faster for simple problems, while the "Conjugate Gradient" guide was better for very complex, tricky mountains.

Summary

In short, this paper introduces a new, smarter way to solve giant 3D math puzzles. By packing the data efficiently (Tucker format), using better hiking guides (Subspace Gradient Descent), and giving them a GPS (Preconditioning) and a fast-forward button (Mixed Precision), the authors have created tools that solve complex scientific problems much faster than the current best methods. This helps scientists simulate everything from quantum chemistry to financial risks more quickly and cheaply.

1. Problem Statement

The paper addresses the numerical solution of linear tensor equations of the form:
$\mathcal{A}(\mathbf{x}) = \mathbf{c}$
where:

$\mathcal{A}$ is a symmetric and positive-definite multilinear operator acting on the tensor space $\mathbb{R}^{n_1 \times \dots \times n_d}$ .
$\mathbf{x}$ and $\mathbf{c}$ are tensors of order $d$ .
The operator $\mathcal{A}$ is often a sum of tensorized operators (e.g., arising from the discretization of 3D partial differential equations with separable coefficients).

Challenges:

Solving these equations directly is computationally prohibitive due to the "curse of dimensionality" (the size of the full tensor grows exponentially with $d$ ).
Existing solvers often lack convergence guarantees (optimization-based methods) or struggle with memory management and parallelization (classical iterative methods).
The Alternating Minimal Energy (AMEn) algorithm is currently a state-of-the-art solver but operates in the Tensor-Train (TT) format, which may be less efficient for problems with a modest number of modes ( $d \approx 3$ ).

2. Methodology

The authors propose two new gradient-descent type methods adapted to the Tucker tensor format, generalizing the Subspace Conjugate Gradient (Ss-cg) method previously developed for matrices.

A. Core Framework: Subspace Gradient Descent

The methods approximate the solution $\mathbf{x}$ iteratively within a subspace. The update rule is:
$\mathbf{x}_{k+1} = \mathbf{x}_k + \mathcal{P}_k \boldsymbol{\alpha}_k$
where $\mathcal{P}_k$ is a tensor operator defining the search subspace, and $\boldsymbol{\alpha}_k$ is a core tensor coefficient determined by minimizing the quadratic energy functional $\Phi(\mathbf{x}) = \frac{1}{2}\langle \mathbf{x}, \mathcal{A}(\mathbf{x}) \rangle - \langle \mathbf{x}, \mathbf{c} \rangle$ within the subspace.

The minimization leads to a reduced multilinear system of much smaller dimensions:
$\mathcal{P}_k^T \mathcal{A} \mathcal{P}_k \boldsymbol{\alpha}_k = \mathcal{P}_k^T \mathbf{r}_k$
where $\mathbf{r}_k$ is the residual.

B. The Two Proposed Algorithms

Tk–ss–cg (Tensor Subspace Conjugate Gradient):
- Generalizes the Conjugate Gradient method.
- The search direction $\mathcal{P}_{k+1}$ is updated using a recurrence involving the previous direction and the new residual, ensuring $\mathcal{A}$ -orthogonality (conjugacy).
- Requires solving for two coefficient tensors ( $\boldsymbol{\alpha}_k$ and $\boldsymbol{\beta}_k$ ) at each step.
Tk–ss–sd (Tensor Subspace Steepest Descent):
- Generalizes the Steepest Descent method.
- The search direction is simply the current residual ( $\mathcal{P}_{k+1} \propto \mathbf{r}_{k+1}$ ).
- Avoids solving for $\boldsymbol{\beta}_k$ , reducing computational cost per iteration, though it may require more iterations to converge.

C. Implementation Strategies

To handle the high dimensionality and maintain efficiency, the authors employ several key techniques:

Tucker Format: Tensors are represented as a core tensor multiplied by factor matrices ( $\mathbf{x} = \mathcal{C} \times_1 U_1 \times_2 \dots \times_d U_d$ ). This is chosen as it is highly competitive for low-mode ( $d=3$ ) problems compared to TT format.
Rank Truncation: Operations (sums and products) increase tensor rank. The authors use ST-HOSVD (Sequentially Truncated Higher-Order SVD) and truncated SVDs to round tensors back to a maximum allowed multilinear rank, controlling memory usage.
Mixed-Precision Computing: To reduce costs, the solution of the reduced systems (for $\boldsymbol{\alpha}_k$ and $\boldsymbol{\beta}_k$ ) and preconditioner applications are performed in single precision, while the main iterative loop remains in double precision.
Preconditioning: Several strategies are implemented to accelerate convergence:
- P–Eig: Uses eigenvalue factorization of the Laplacian operator.
- P–FFT: Uses Fast Fourier Transforms (specifically Discrete Sine Transforms) to exploit the structure of the Laplacian, offering very fast application.
- P–InnOut: An inner-outer loop where a few iterations of the proposed method solve the preconditioning system.

3. Key Contributions

Generalization of Ss-cg: Successfully extends the subspace gradient descent framework from matrices to high-order tensors using the Tucker format.
Algorithmic Variants: Introduces both a Conjugate Gradient variant (Tk–ss–cg) and a Steepest Descent variant (Tk–ss–sd), offering a trade-off between per-iteration cost and convergence speed.
Efficient Implementation: Develops specific routines for Tucker-rounding-sum and Tucker-rounding-product to manage rank growth without exploding memory.
Preconditioning for Tensors: Proposes and tests effective preconditioners (P–FFT, P–Eig) tailored for tensor operators arising from PDEs, demonstrating that P–FFT is particularly effective for Laplacian-like problems.
Mixed-Precision Strategy: Demonstrates that solving reduced sub-problems in single precision significantly reduces computational time without sacrificing final accuracy.

4. Experimental Results

The methods were tested on 3D discretizations of elliptic PDEs (Poisson equation, diffusion with variable coefficients, reaction-diffusion, and discontinuous coefficients) and compared against the AMEn algorithm (which uses the TT format).

Performance vs. AMEn:
- Unpreconditioned: The proposed methods are competitive with AMEn. Tk–ss–sd often has lower CPU time for moderate tolerances due to lower per-iteration cost, despite requiring slightly more iterations.
- Preconditioned: The preconditioned versions (especially P–FFT Tk–ss–sd) significantly outperform AMEn.
  - In Example 5.1 (Poisson), P–FFT Tk–ss–sd solved the problem in 2–3 iterations with CPU times of 0.05–0.06 seconds, compared to AMEn's 0.23–0.75 seconds.
  - The P–FFT preconditioner is highly effective because the problem operators are spectrally close to the Laplacian.
Scalability: The methods scale well with grid size ( $n=500$ to $n=1000$ ).
Rank Sensitivity: The convergence of Tk–ss–sd slows down as iterations increase if the maximum rank is too low, but it tracks Tk–ss–cg closely when sufficient rank is allowed.
Precision: Mixed-precision strategies successfully reduced costs. P–FFT required double precision in some cases to avoid stagnation, but the overhead was negligible due to the extremely low iteration count.

5. Significance

Competitiveness: The paper establishes that subspace gradient descent methods in the Tucker format are a robust and often superior alternative to the dominant AMEn algorithm for 3D tensor problems (modest $d$ ).
Efficiency: By combining the Tucker format (efficient for low $d$ ) with mixed-precision arithmetic and specialized preconditioners, the authors achieve significant speedups.
Practicality: The approach offers a transparent management of memory and accuracy, making it suitable for large-scale scientific computing applications involving 3D PDEs, quantum chemistry, and imaging where tensor equations naturally arise.
Future Direction: The work highlights the potential of subspace projection methods for tensors and suggests that further tuning of preconditioners could yield even greater performance gains.