Eigen-Contrib: Eigen::BartelsStewart< KroneckerSumType

Searching...

No Matches

#include <contrib/Eigen/src/StructuredMatrices/KroneckerSum.h>

Detailed Description

template<typename KroneckerSumType>
class Eigen::BartelsStewart< KroneckerSumType >

Direct solver for Kronecker-sum systems \( (A_1 \oplus A_2 \oplus \cdots \oplus A_d)\,x = b \).

Flattening the (possibly nested) KroneckerSum into its non-sum factors \( A_k = Q_k T_k Q_k^H \) (complex Schur forms, \( Q_k \) unitary, \( T_k \) upper triangular) gives

\[ A_1 \oplus \cdots \oplus A_d = Q\,(T_1 \oplus \cdots \oplus T_d)\,Q^H, \qquad Q = Q_1 \otimes \cdots \otimes Q_d, \]

so a solve applies \( Q^H \), back-substitutes the upper triangular \( T_1 \oplus \cdots \oplus T_d \) and applies \( Q \). Both transforms are one product per factor along its own index; the back substitution recurses over the factors with accumulated shifts,

\[ \big((\sigma + T_1(i,i)) I + T_2 \oplus \cdots \oplus T_d\big)\,y_i = c_i - \textstyle\sum_{j > i} T_1(i,j)\,y_j, \]

the Bartels-Stewart algorithm [1] for two factors. When every factor is exactly Hermitian, the Schur forms are the eigendecompositions, \( T_k \) is real diagonal, and the triangular solve is a division by the eigenvalue sums – the fast diagonalization method [2], in real arithmetic for real factors. Real factors that are not all symmetric take real Schur forms instead, for up to four factors: \( Q_k \) is orthogonal and \( T_k \) quasi-upper-triangular, with a 2x2 diagonal block per complex conjugate eigenvalue pair. The block rows of the recurrence are then solved jointly: with \( \Sigma \) the Kronecker sum of the diagonal blocks chosen at the outer factors (at most \( 2^{d-1} \) wide) in place of \( \sigma \),

\[ \big((\Sigma \oplus T_1(I,I)) \oplus T_2 \oplus \cdots \oplus T_d\big)\,Y_I = C_I - \textstyle\sum_{J > I} Y_J\,T_1(I,J)^T \]

for each diagonal block \( I \) of \( T_1 \), ending in dense solves at most \( 2^d \) wide, as in LAPACK's xTRSYL for two factors. More real factors, and complex ones, take the complex Schur form. Either way the setup costs one \( O(n_k^3) \) decomposition per factor, each solve \( O(N \sum_k n_k) \) per right-hand side, \( N = \prod_k n_k \), plus up to \( O(4^d N) \) for the coupled solves of the real path; sparse factors are densified for the decomposition. transpose().solve() and adjoint().solve() reuse the decompositions:

\[ (A_1 \oplus \cdots \oplus A_d)^H = Q\,(T_1^H \oplus \cdots \oplus T_d^H)\,Q^H \]

is solved with the same transforms around a forward substitution (on the complex path, on copies of the \( T_k^H \) that compute keeps), and \( M^T x = b \) as \( M^H \bar x = \bar b \).

The system is singular exactly when some sum \( \lambda_{i_1}(A_1) + \cdots + \lambda_{i_d}(A_d) \) vanishes. As with PartialPivLU nothing detects it: the computed sum is typically of order \( \epsilon \sum_k \|A_k\| \) rather than zero, and the solution huge; it is non-finite when a sum vanishes exactly, as it can for diagonal, triangular or identity factors, whose decompositions are exact.

The solve is backward stable relative to the factors, with residual \( \|b - Mx\| = O\big((\sum_k n_k)\,\epsilon\,(\sum_k \|A_k\|)\,\|x\|\big) \). That is relative to \( \sum_k \|A_k\| \), not \( \|M\| \); with a shift split across the factors, as in \( (A + cI) \oplus (B - cI) \), the first grows with \( c \) and the second does not.

info() reports InvalidInput for a non-finite factor and NoConvergence when a Schur or eigenvalue iteration fails; the solve then returns NaN.

Template Parameters

KroneckerSumType the KroneckerSum type to solve with.

See also: class KroneckerSum

Inheritance diagram for Eigen::BartelsStewart< KroneckerSumType >: