 |
Eigen-Contrib
5.0.1
|
|
|
Eigen-Contrib
5.0.1
|
|
#include <contrib/Eigen/src/StructuredMatrices/KroneckerOperator.h>
The Kronecker product \( A \otimes B \) as an implicit operator that is never materialized.
For A of size m1 x n1 and B of size m2 x n2, the Kronecker product is the m1*m2 x n1*n2 block matrix whose block (i,j) is A(i,j)*B.
Throughout, \( \mathrm{vec}(X) \) stacks the columns of X, whatever the storage order of the operands and independently of EIGEN_DEFAULT_TO_ROW_MAJOR: for X of size n2 x n1, entry j*n2+i of \( \mathrm{vec}(X) \) is X(i,j), and its inverse \( \mathrm{mat}(x) \) is x.reshaped(n2,n1). Under this convention the product is \( (A \otimes B)\,\mathrm{vec}(X) = \mathrm{vec}(B X A^T) \); with row-major stacking \( \mathrm{vec}_r \) of Z of size n1 x n2 the same product reads \( (A \otimes B)\,\mathrm{vec}_r(Z) = \mathrm{vec}_r(A Z B^T) \).
This class stores only the two factors and evaluates every operation through them:
A and B ( \( (A \otimes B)^{-1} = A^{-1} \otimes B^{-1} \)); minimum-norm least-squares solves ( \( (A \otimes B)^+ = A^+ \otimes B^+ \)) use one complete orthogonal decomposition per factor, deciding each factor's rank on its own; rank goes through the factor SVDs, thresholding the pairwise singular-value products \( \sigma_i(A)\,\sigma_j(B) \) – the singular values of the Kronecker product – at the product level;KroneckerOperator objects themselves, never materialized;The class is closed under transpose, conjugate and adjoint ( \( (A \otimes B)^T = A^T \otimes B^T \)). operator* returns an Eigen product expression, so the operator plugs into the matrix-free iterative solvers, and it can be assigned to a dense matrix when an explicit representation is needed. As with any matrix-free operator, the iterative solvers must be instantiated with IdentityPreconditioner (e.g. ConjugateGradient<KroneckerOperator<MatrixXd,MatrixXd>,Lower|Upper,IdentityPreconditioner>): the default preconditioners read individual coefficients through col() or InnerIterator, which the structured operators do not expose.
In contrast to kroneckerProduct() (the KroneckerProduct module), which builds an expression meant to be evaluated into a dense matrix, this class is an operator meant to be applied and solved with, without ever forming the product.
Either factor may be a DiagonalMatrix. A diagonal factor is stored as its diagonal – O(n) instead of O(n^2) – its side of every product is a diagonal scaling instead of a GEMM, solve divides entrywise instead of factorizing, and transpose, conjugate, adjoint, inverse and determinant never leave diagonal form.
Either factor may also be an identity, passed as an Identity() expression (MatrixXd::Identity(p, p)): makeKroneckerOperator() stores it as an internal identity factor holding only its dimensions (assignable, unlike the expression), its side of a product is skipped, and its solves are the identity map. This covers the identity-Kronecker operators \( I \otimes A \) and \( A \otimes I \) ubiquitous in finite-difference and Sylvester/Lyapunov settings, e.g.
A rectangular Identity(m, n) is the m x n matrix with ones on the main diagonal. The unit DiagonalMatrix VectorXd::Ones(p).asDiagonal() describes the same operator but is applied as a scaling.
The decomposition family (eigenvalues, eigenvectors, singularValues, matrixU, matrixV, leastSquaresSolve, rank) currently materializes a diagonal or identity factor densely for the factor decomposition.
Either factor may also be a SparseMatrix, stored compressed. Its side of a product is a sparse-dense product – O(nnz) instead of O(m n) per column of the reshaped right-hand side – solve factorizes it once with SparseLU and back-substitutes, transpose, conjugate and adjoint stay sparse, and determinant accumulates the SparseLU pivots in the same balanced form as the dense LU path, scaling small factors up before factorization to keep the elimination out of the subnormal range. Its inverse is dense (one SparseLU solve against the identity), and the decomposition family densifies it like a diagonal factor. A product with non-finite data propagates Inf/NaN as the two factor products do, not as a product with the materialized sparse matrix would. A sparse factor whose SparseLU factorization fails – an exactly zero or non-finite pivot column – solves and inverts to NaN, and has determinant 0 when exactly singular (NaN when non-finite). Assigning the operator to a SparseMatrix materializes the product sparsely, every inner vector reserved to its exact size. For a sparse A of size n with nnz(A) stored entries:
Finally, either factor may be a KroneckerOperator itself, so products of three or more factors stay implicit: \( A \otimes (B \otimes C) \) applies \( B \otimes C \) through its own vec identity, and solves, the determinant, the transposition family, materialization, eigenvalues, eigenvectors and the SVD (singularValues, matrixU, matrixV, whose vector matrices nest the same way) recurse into the nested factors; rank and leastSquaresSolve materialize a nested factor. makeKroneckerOperator(a, b, c, ...) builds this right-nested form; e.g. the middle term \( I_p \otimes A \otimes I_q \) of a 3-D finite-difference operator is
| LhsMatrix | the type of the left factor A: a dense Matrix, a DiagonalMatrix to exploit diagonal structure, a SparseMatrix to exploit sparsity, an identity factor (what makeKroneckerOperator() stores an Identity() expression as), a KroneckerOperator, or a KroneckerSum. |
| RhsMatrix | the type of the right factor B, under the same convention; its scalar type must match that of LhsMatrix. |
Inheritance diagram for Eigen::KroneckerOperator< LhsMatrix, RhsMatrix >:Public Member Functions | |
| KroneckerOperator< typename LhsOps::TransposedFactor, typename RhsOps::TransposedFactor > | adjoint () const |
| Scalar | coeff (Index row, Index col) const |
| KroneckerOperator | conjugate () const |
| Scalar | determinant () const |
| ComplexVector | eigenvalues () const |
| KroneckerOperator< typename LhsSpectrum::Eigenvectors, typename RhsSpectrum::Eigenvectors > | eigenvectors () const |
| KroneckerOperator< typename LhsOps::InverseFactor, typename RhsOps::InverseFactor > | inverse () const |
| template<typename LhsDerived, typename RhsDerived> | |
| KroneckerOperator (const EigenBase< LhsDerived > &a, const EigenBase< RhsDerived > &b) | |
| template<typename Rhs> | |
| Matrix< Scalar, ColsAtCompileTime, Rhs::ColsAtCompileTime > | leastSquaresSolve (const MatrixBase< Rhs > &b) const |
| const LhsMatrix & | lhs () const |
| KroneckerOperator< typename LhsSpectrum::SingularVectors, typename RhsSpectrum::SingularVectors > | matrixU () const |
| KroneckerOperator< typename LhsSpectrum::SingularVectors, typename RhsSpectrum::SingularVectors > | matrixV () const |
| template<typename Rhs> | |
| Product< KroneckerOperator, Rhs > | operator* (const MatrixBase< Rhs > &x) const |
| Index | rank () const |
| const RhsMatrix & | rhs () const |
| RealVector | singularValues () const |
| template<typename Rhs> | |
| Matrix< Scalar, ColsAtCompileTime, Rhs::ColsAtCompileTime > | solve (const MatrixBase< Rhs > &b) const |
| KroneckerOperator< typename LhsOps::TransposedFactor, typename RhsOps::TransposedFactor > | transpose () const |
|
inline |
Builds the operator A (x) B from the two factors, evaluated into the operator's factor types: a dense expression into a Matrix, a diagonal one into a DiagonalMatrix (stored as its diagonal), a sparse one into a compressed SparseMatrix, an Identity() expression into an identity factor holding its dimensions; a KroneckerOperator is copied as it is.
|
inline |
|
inline |
|
inline |
|
inline |
A of size n1 and B of size n2. The product is accumulated from the factor LU diagonals (the SparseLU pivots for a sparse factor, the diagonal itself for a diagonal factor, skipping the LU; 1 for an identity; recursively for a nested Kronecker factor, whose own factors must then be square; the LU of the materialized matrix for a KroneckerSum factor) in the balanced form m * 2^e – every factor and the running product are renormalized to unit magnitude with the power of two tracked separately – so the partial products (in particular det(A) and det(B) themselves, which can overflow or underflow on their own) never leave the representable range when the determinant itself is representable.
|
inline |
i*n2 + j is \( \lambda_i(A)\,\mu_j(B) \), matching column i*n2 + j of eigenvectors. The set is not sorted – there is no canonical eigenvalue order, and sorting would break the Kronecker structure of the eigenvector matrix. A nested Kronecker or KroneckerSum factor contributes the eigenvalues of its own factors – the factors of every nested Kronecker product must then be square – with the accuracy KroneckerSum::eigenvalues() describes for a sum.
|
inline |
i*n2 + j is \( v_i(A) \otimes v_j(B) \) and matches eigenvalues()[i*n2 + j]. Assign it to a dense matrix to materialize.
|
inline |
SparseLU solve against the identity (NaN when the factorization fails).
|
inline |
The numerical rank is decided per factor, from the diagonal of the pivoted QR inside CompleteOrthogonalDecomposition (so, as for any column-pivoted QR, a near-deficiency of the Kahan type can go undetected). This is backward stable: the decompositions perturb A and B separately, so each product singular value \( \sigma_i(A)\,\sigma_j(B) \) inherits only its factors' relative errors. It therefore need not agree with rank, which thresholds the pairwise singular-value products at the product level: factors that are each full rank can form modes that rank drops but this method inverts.
Each factor is scaled by an exact power of two before its decomposition, so factor magnitudes alone cannot over- or underflow the decompositions or the intermediates; the right-hand side is not rescaled. A non-finite factor solves to NaN.
|
inline |
A.
|
inline |
i*k_B + j matches singularValues()[i*k_B + j].
|
inline |
i*k_B + j matches singularValues()[i*k_B + j].
|
inline |
|
inline |
min(rows(),cols()) * epsilon * sigma_max(A) * sigma_max(B) (the SVDBase convention), and that reach the smallest normal number (the SVDBase threshold clamp, so subnormal products count as exact zeros). The relative comparison is made in ratio space, \( (\sigma_i(A)/\sigma_{max}(A))(\sigma_j(B)/\sigma_{max}(B)) \) against min(rows(),cols()) * epsilon, and the clamp in exponent space, so that neither the thresholds nor the products can spuriously under- or overflow. Thresholding the products matters: factors that are each full rank against their own threshold can still form pairwise products that are negligible at the product level, so the rank can be smaller than the product of the factor ranks. This costs two SVDs; leastSquaresSolve decides each factor's rank on its own instead.
|
inline |
B.
|
inline |
i*k_B + j is \( \sigma_i(A)\,\sigma_j(B) \) with k_A, k_B the factor thin ranks min(rows, cols), matching the columns of matrixU and matrixV. The values are not sorted (sorting would break the Kronecker structure of U and V); for rectangular shapes the full SVD pads this set with min(rows(),cols()) - k_A*k_B structural zeros.
|
inline |
SparseLU per sparse factor (a diagonal factor is solved by entrywise division instead): reshaping b column-wise as mat(b) of size n2 x n1, the system reads \( B X A^T = \mathrm{mat}(b) \), so \( X = B^{-1} \mathrm{mat}(b) A^{-T} \). Right-hand sides are solved in cache-sized batches, one solve per factor per batch, at O(n1^3 + n2^3 + nrhs (n1 + n2) n1 n2) total cost for dense factors; a diagonal factor contributes only O(nrhs n1 n2), a sparse one its factorization plus its substitutions. As in any chain of solves, \( B^{-1} \mathrm{mat}(b) \) can overflow (B tiny) or underflow (B huge, A tiny) on extreme factor magnitudes even when the solution is representable. PartialPivLU: a singular dense factor substitutes Inf/NaN through the solution, and a sparse factor whose SparseLU factorization fails solves to NaN. Use leastSquaresSolve for rank-deficient or rectangular factors.
|
inline |
1.13.2