Multi-Subspace Power Method (MSPM)

• MSPM is a tensor decomposition algorithm that recovers low-rank, partially symmetric CP decompositions from real tensors with arbitrary partial symmetry.
• It employs an orthonormal-slice transformation combined with a shifted higher-order power iteration to extract partially symmetric singular-vector tuples.
• MSPM demonstrates superior accuracy and runtime performance compared to methods like ALS and nonlinear least-squares, with guarantees for global and local linear convergence.

The Multi-Subspace Power Method (MSPM) is a tensor decomposition algorithm for recovering low-rank, partially symmetric CP decompositions from real tensors exhibiting arbitrary partial symmetry—ranging from fully symmetric to fully asymmetric cases. MSPM generalizes the subspace power method beyond symmetric tensors by introducing an orthonormal-slice transformation and a globally convergent, shifted higher-order power iteration for extracting partially symmetric singular-vector tuples (pSVTs). The method facilitates efficient and accurate recovery of tensor factors and demonstrates theoretical and empirical advantages over classical approaches(Wang et al., 21 Oct 2025).

1. Problem Formulation and Preliminaries

Consider positive integers $d_1,\dots,d_\ell$, $m_1,\dots,m_\ell$, and $r$. Let the input tensor be

$$\mathcal{T} \in S^{d_1}(\mathbb{R}^{m_1}) \otimes \cdots \otimes S^{d_\ell}(\mathbb{R}^{m_\ell}),$$

where $S^{d_i}(\mathbb{R}^{m_i})$ denotes the space of order-$d_i$ symmetric tensors over $\mathbb{R}^{m_i}$. The target is the partially symmetric CP decomposition: \begin{equation} \mathcal{T} = \sum_{i=1}^{r} \lambda_i\, (v_i^{(1)})^{\otimes d_1} \otimes \cdots \otimes (v_i^{(\ell)})^{\otimes d_\ell}, \end{equation} with vectors $v_i^{(j)} \in \mathbb{S}^{m_j-1}$ and $\lambda_i \in \mathbb{R}\setminus\{0\}$. The decomposition problem is generally NP-hard, but for low $r$ and generic tensors, algorithmic solutions exist that can guarantee recovery(Wang et al., 21 Oct 2025).

2. Orthonormal-Slice Transformation

MSPM starts by choosing a flattening pattern $m_1,\dots,m_\ell$0, $m_1,\dots,m_\ell$1, and computes the $m_1,\dots,m_\ell$2-flattening of $m_1,\dots,m_\ell$3: $m_1,\dots,m_\ell$4 with

$m_1,\dots,m_\ell$5

A thin singular value decomposition (SVD) yields leading $m_1,\dots,m_\ell$6 left singular vectors, which are reshaped into a tensor $m_1,\dots,m_\ell$7 of dimension

$m_1,\dots,m_\ell$8

with orthonormal slices along the last ($m_1,\dots,m_\ell$9th) mode. The original rank-one tensor summands map to partially symmetric rank-one components in $r$0 with shared factors $r$1 and a standard basis vector in $r$2(Wang et al., 21 Oct 2025).

3. Partially Symmetric Singular-Vector Tuples

A partially symmetric singular-vector tuple (pSVT) is a tuple $r$3 and scalar $r$4 solving, for each $r$5, \begin{equation} \mathcal{A} \cdot (v^{(1)})^{\otimes(a_1-\delta_{1j})} \otimes\cdots\otimes (v^{(\ell)})^{\otimes(a_\ell-\delta_{\ell j})} = \sigma\, v^{(j)}, \end{equation} where $r$6 is a partially symmetric tensor and $r$7 is the Kronecker delta. For the orthonormal-slice tensor $r$8, all pSVTs satisfy $r$9; those with $$\mathcal{T} \in S^{d_1}(\mathbb{R}^{m_1}) \otimes \cdots \otimes S^{d_\ell}(\mathbb{R}^{m_\ell}),$$0 correspond exactly to the rank-one components of the decomposition(Wang et al., 21 Oct 2025).

4. The Multi-Subspace Power Method: Algorithmic Structure

The MSPM algorithm proceeds as follows:

1. Subspace Extraction: Compute the $$\mathcal{T} \in S^{d_1}(\mathbb{R}^{m_1}) \otimes \cdots \otimes S^{d_\ell}(\mathbb{R}^{m_\ell}),$$1-flattening $$\mathcal{T} \in S^{d_1}(\mathbb{R}^{m_1}) \otimes \cdots \otimes S^{d_\ell}(\mathbb{R}^{m_\ell}),$$2 of $$\mathcal{T} \in S^{d_1}(\mathbb{R}^{m_1}) \otimes \cdots \otimes S^{d_\ell}(\mathbb{R}^{m_\ell}),$$3, perform SVD and truncate to top $$\mathcal{T} \in S^{d_1}(\mathbb{R}^{m_1}) \otimes \cdots \otimes S^{d_\ell}(\mathbb{R}^{m_\ell}),$$4 singular vectors; construct $$\mathcal{T} \in S^{d_1}(\mathbb{R}^{m_1}) \otimes \cdots \otimes S^{d_\ell}(\mathbb{R}^{m_\ell}),$$5 with orthonormal slices.
2. Rank-one Recovery: For $$\mathcal{T} \in S^{d_1}(\mathbb{R}^{m_1}) \otimes \cdots \otimes S^{d_\ell}(\mathbb{R}^{m_\ell}),$$6 to $$\mathcal{T} \in S^{d_1}(\mathbb{R}^{m_1}) \otimes \cdots \otimes S^{d_\ell}(\mathbb{R}^{m_\ell}),$$7, apply the shifted higher-order power method (PS-HOPM) to $$\mathcal{T} \in S^{d_1}(\mathbb{R}^{m_1}) \otimes \cdots \otimes S^{d_\ell}(\mathbb{R}^{m_\ell}),$$8 to find a pSVT $$\mathcal{T} \in S^{d_1}(\mathbb{R}^{m_1}) \otimes \cdots \otimes S^{d_\ell}(\mathbb{R}^{m_\ell}),$$9.
3. Completion: Extend recovered vectors to the full tensor rank-one term.
4. Coefficient Recovery: Calculate the $S^{d_i}(\mathbb{R}^{m_i})$0 via projection.
5. Deflation and Update: Remove the recovered component via a one-rank SVD downdate.
6. Termination: Repeat until $S^{d_i}(\mathbb{R}^{m_i})$1 components are recovered.

The primary subroutine, PS-HOPM, performs alternating block updates with shift parameters $S^{d_i}(\mathbb{R}^{m_i})$2 to ensure convexity, followed by normalization, and updates the $S^{d_i}(\mathbb{R}^{m_i})$3 value as a fully contracted tensor product(Wang et al., 21 Oct 2025).

5. Convergence Properties

• Global Convergence: Given shift parameters $S^{d_i}(\mathbb{R}^{m_i})$4, the shifted objective,

$S^{d_i}(\mathbb{R}^{m_i})$5

is maximized by alternating block ascent, converging to pSVTs from any initialization(Wang et al., 21 Oct 2025).

• Local Linear Convergence: In the special case $S^{d_i}(\mathbb{R}^{m_i})$6, every strict local maximizer is a fixed point for PS-HOPM under which convergence is locally linear. This is established via spectral analysis of the Jacobian at the fixed point.

6. Computational Complexity and Empirical Results

Let $S^{d_i}(\mathbb{R}^{m_i})$7 be the total number of tensor entries.

• Orthonormal-Slice Construction: The thin SVD step has $S^{d_i}(\mathbb{R}^{m_i})$8 cost.
• PS-HOPM Iteration: Each update costs $S^{d_i}(\mathbb{R}^{m_i})$9 arithmetic, relying on two high-dimensional contractions.
• Deflation: SVD downdate per step requires $d_i$0.
• Total Cost: $d_i$1.

Empirically, for tensors up to size $d_i$2 and rank $d_i$3, MSPM achieves higher accuracy and runtime performance than Alternating Least Squares (ALS), nonlinear least-squares solvers (TensorLab NLS/MINF), simultaneous diagonalization (FFDIAG, Jacobi, Jennrich), and generalized SVD (HOGSVD). In log-error versus log-time plots, MSPM occupies the best-accuracy, best-time region(Wang et al., 21 Oct 2025).

7. Limitations and Directions for Generalization

• Rank Constraints: The theoretical recovery guarantee holds when $d_i$4; above this, uniqueness is lost and performance may degrade.
• Exploitation of Symmetry: MSPM leverages tensor partial symmetry via flattenings tailored to the symmetry type. In the fully asymmetric case ($d_i$5), it reduces to a deflationary CP power method, still performing competitively in experiments.
• Unresolved Issues: Open questions include extending local linear convergence to $d_i$6 blocks, more precise robustness/nonideal SVD truncation analysis, and adaptive shift strategies for accelerated PS-HOPM convergence without tight spectral norm estimation.

A plausible implication is that MSPM unifies and subsumes a range of tensor decomposition strategies under a common algebraic-geometric framework, facilitating future advances in the study of low-rank tensor varieties and algorithms(Wang et al., 21 Oct 2025).