Factorization methods

Many forms of matrix factorization such as eigenvalue decomposition, LU factorization, Cholesky factorization etc., can be computed efficiently. The decomposition of the Kronecker product is the Kronecker product of the decompositions. We have overloaded some of the factorization functions from LinearAlgebra to compute the factorization of instances of AbstractKroneckerProduct.

Eigenvalue decomposition

The function eigen of LinearAlgebra is overloaded to compute the decomposition of AbstractKroneckerProducts. The result is a factorization of the Eigen type, containing a vector of the eigenvalues and a matrix with the eigenvectors. Just like long-time users would expect! The eigenvectors are structured as Kronecker products and can be processed accordingly.

The functions det, logdet, inv and \ are overloaded the make use of this decomposition.

The eigenvalue decomposition of matrices can be used to solve large systems of the form:

\[(A \otimes B + c\cdot I) \mathbf{x} = \mathbf{b}\]

The case where $A$ and $B$ are positive semi-definite frequently occurs in machine learning, for example in ridge regression.


julia> A, B = rand(10, 10), randn(4, 4);
julia> As, Bs = (A, B) .|> X -> X * X'; # make positive definite([2.923612267419591 1.458953187121179 … 2.6134296301570794 2.921623438268228; 1.458953187121179 1.2348238914394167 … 1.6932109656424656 1.3966044795930572; … ; 2.6134296301570794 1.6932109656424656 … 4.012084917238429 2.9485603660442985; 2.921623438268228 1.3966044795930572 … 2.9485603660442985 4.138589040805155], [8.212216560569244 0.35959547703244815 3.784715273838838 -6.878503127866228; 0.35959547703244815 5.800866391145329 0.13313135643791735 -2.4334768227778314; 3.784715273838838 0.13313135643791735 4.376486369440425 1.7168499108973525; -6.878503127866228 -2.4334768227778314 1.7168499108973525 18.52088887897743])
julia> K = As ⊗ Bs40×40 Kronecker.KroneckerProduct{Float64, Matrix{Float64}, Matrix{Float64}}: 24.0093 1.05132 11.065 … 1.0506 11.0575 -20.0964 1.05132 16.9595 0.389224 16.9479 0.38896 -7.1097 11.065 0.389224 12.7951 0.38896 12.7864 5.01599 -20.1101 -7.11454 5.0194 -7.1097 5.01599 54.1111 11.9812 0.524633 5.52172 0.502213 5.28575 -9.60655 0.524633 8.46319 0.194232 … 8.10152 0.185932 -3.3986 5.52172 0.194232 6.38509 0.185932 6.11222 2.39776 -10.0354 -3.55033 2.5048 -3.3986 2.39776 25.8664 17.3236 0.758562 7.98381 0.968094 10.1891 -18.5181 0.758562 12.2369 0.280839 15.6169 0.358413 -6.55135 ⋮ ⋱ -18.8725 -6.6767 4.7105 -7.86884 5.55157 59.8887 21.4621 0.939777 9.89109 1.06029 11.1595 -20.2817 0.939777 15.1602 0.347929 17.1042 0.392546 -7.17525 9.89109 0.347929 11.4376 0.392546 12.9043 5.06224 -17.9765 -6.35972 4.48687 … -7.17525 5.06224 54.61 23.993 1.0506 11.0575 1.48822 15.6634 -28.4673 1.0506 16.9479 0.38896 24.0074 0.550976 -10.0712 11.0575 0.38896 12.7864 0.550976 18.1125 7.10534 -20.0964 -7.1097 5.01599 -10.0712 7.10534 76.6503
julia> E = eigen(K)Eigen{Float64, Float64, Kronecker.KroneckerProduct{Float64, Matrix{Float64}, Matrix{Float64}}, Vector{Float64}} values: 40-element Vector{Float64}: 9.060463112706263e-5 0.001445827914292918 0.0023247809976477698 0.005890155378432289 0.010013955507711825 0.15979819381674532 0.2569433061646792 0.6510015344625429 0.05636829423551302 0.8994998630090828 ⋮ 24.98689239884551 0.7171049716952264 11.44323830534869 18.39985429859397 46.618585092054424 8.806456677500273 140.52947108724626 225.96066984423558 572.5027244045117 vectors: 40×40 Kronecker.KroneckerProduct{Float64, Matrix{Float64}, Matrix{Float64}}: -0.0546617 -0.00376689 0.0582939 … 0.204659 0.137353 -0.0104751 0.086632 -0.0126244 -0.0443221 0.0439426 0.0635929 0.0163425 0.0601591 0.211208 0.0027591 -0.0280924 0.0120207 0.0274625 0.0964159 -0.277399 0.476591 0.0328432 -0.508259 0.126027 0.0845803 0.0913315 -0.755337 0.110071 … -0.027293 0.0270593 -0.554461 -0.142489 -0.524522 0.130059 0.00169902 0.244935 -0.104808 -0.239443 0.0593717 -0.170819 -0.0794148 -0.0054727 0.0846918 0.168175 0.112867 -0.0152187 0.125863 -0.0183413 -0.0364209 0.0361091 ⋮ ⋱ -0.0607493 0.0259946 0.0593873 0.114001 -0.327993 0.100716 0.00694063 -0.107409 0.232861 0.15628 0.0193007 -0.159622 0.023261 -0.0504296 0.0499978 -0.117172 -0.0301116 -0.110845 0.240312 0.0031393 0.0517612 -0.0221486 -0.0506007 … 0.109702 -0.315624 0.0548171 0.0037776 -0.0584596 0.232463 0.156013 0.0105049 -0.0868782 0.0126603 -0.0503435 0.0499124 -0.0637736 -0.0163889 -0.0603301 0.239901 0.00313394 0.0281722 -0.0120549 -0.0275406 0.109514 -0.315085
julia> logdet(E)10.354639610069293
julia> b = randn(40);
julia> (E + 0.1I) \ b # solve a system40-element Vector{Float64}: 5.026081572604652 -0.9803398975448284 -7.549796327360351 3.523446931054747 7.917517240208104 -9.118733947227131 3.0624949888104696 -4.681203682086042 -1.538121710532522 3.9403292243582952 ⋮ 1.4901421800748147 5.147587308714193 -0.4491420967583311 -3.5969240171074235 1.1158912544047401 -1.2440461139619112 -2.9040372634471576 2.0954873641221696 -1.5254509168222792
LinearAlgebra.eigenFunction
eigen(K::AbstractKroneckerProduct)

Wrapper around eigen from the LinearAlgebra package. If the matrices of an AbstractKroneckerProduct instance are square, performs Eigenvalue decompositon on them and returns an Eigen type. Otherwise, it collects the instance and runs eigen on the full matrix. The functions, \, inv, and logdet are overloaded to efficiently work with this type.

source
Missing docstring.

Missing docstring for +(E::Eigen, B::UniformScaling). Check Documenter's build log for details.

Missing docstring.

Missing docstring for +(::Eigen, ::UniformScaling). Check Documenter's build log for details.

LinearAlgebra.logdetMethod
logdet(K::Eigen)

Compute the logarithm of the determinant of the eigenvalue decomp of a Kronecker product.

source
Base.invMethod
inv(K::Eigen)

Compute the inverse of the eigenvalue decomp of a Kronecker product. Returns another type of Eigen.

source

Cholesky factorization

Similar to the eigenvalue decomposition, cholesky has been overloaded to allow for efficient Cholesky decomposition of Kronecker products of symmetric and positive definite matrices.


julia> A, B = rand(10, 10), randn(4, 4);
julia> As, Bs = (A, B) .|> X -> X * X'; # make positive definite([3.6506575737985525 2.4389563953188054 … 2.9200058682157524 3.1728907118055143; 2.4389563953188054 2.462263854960425 … 2.012835187347676 2.6123594700089896; … ; 2.9200058682157524 2.012835187347676 … 4.105820790369509 2.8452806642289366; 3.1728907118055143 2.6123594700089896 … 2.8452806642289366 3.8698554581922293], [19.27596400487266 0.8861106032061121 -2.7965035664541404 2.5158843913004345; 0.8861106032061121 4.502179296573766 0.8690381469249804 0.049861905209955246; -2.7965035664541404 0.8690381469249804 3.3102467493373484 0.5098939982244883; 2.5158843913004345 0.049861905209955246 0.5098939982244883 0.7257550919553055])
julia> K = As ⊗ Bs40×40 Kronecker.KroneckerProduct{Float64, Matrix{Float64}, Matrix{Float64}}: 70.3699 3.23489 -10.2091 … 2.81153 -8.873 7.98263 3.23489 16.4359 3.17256 14.2849 2.75736 0.158206 -10.2091 3.17256 12.0846 2.75736 10.5031 1.61784 9.18463 0.182029 1.86145 0.158206 1.61784 2.30274 47.0132 2.16119 -6.82055 2.31484 -7.30547 6.57239 2.16119 10.9806 2.11955 … 11.7613 2.27024 0.130257 -6.82055 2.11955 8.07355 2.27024 8.64755 1.33203 6.13613 0.121611 1.24361 0.130257 1.33203 1.89593 54.5507 2.50768 -7.91406 2.31648 -7.31065 6.57705 2.50768 12.7411 2.45936 11.7696 2.27185 0.13035 ⋮ ⋱ 4.82631 0.0956518 0.978147 0.119003 1.21694 1.73212 56.2859 2.58745 -8.16581 2.52123 -7.95684 7.1584 2.58745 13.1464 2.5376 12.81 2.47266 0.141871 -8.16581 2.5376 9.66594 2.47266 9.41858 1.45079 7.3464 0.145597 1.48889 … 0.141871 1.45079 2.06498 61.1605 2.81153 -8.873 3.42912 -10.8221 9.73611 2.81153 14.2849 2.75736 17.4228 3.36305 0.192958 -8.873 2.75736 10.5031 3.36305 12.8102 1.97322 7.98263 0.158206 1.61784 0.192958 1.97322 2.80857
julia> C = cholesky(K)40×40 CholeskyKronecker{Cholesky{Float64, Matrix{Float64}}, Cholesky{Float64, Matrix{Float64}}} U factor: 40×40 Kronecker.KroneckerProduct{Float64, UpperTriangular{Float64, Matrix{Float64}}, UpperTriangular{Float64, Matrix{Float64}}}: 8.38868 0.385625 -1.21701 … 0.335158 -1.05774 0.951595 0.0 4.03574 0.902404 3.50758 0.784305 -0.051726 0.0 0.0 3.12876 0.0 2.71929 0.902151 0.0 0.0 0.0 0.0 0.0 0.528483 0.0 0.0 -0.0 0.108941 -0.343809 0.309309 0.0 0.0 0.0 … 1.14011 0.254933 -0.0168132 0.0 0.0 0.0 0.0 0.883887 0.293238 0.0 0.0 0.0 0.0 0.0 0.171779 0.0 0.0 -0.0 0.0179664 -0.0567008 0.051011 0.0 0.0 0.0 0.188027 0.0420433 -0.00277281 ⋮ ⋱ 0.0 0.0 0.0 -0.0 -0.0 -0.0151971 0.0 0.0 -0.0 -0.141624 0.446955 -0.402105 0.0 0.0 0.0 -1.48216 -0.331415 0.0218573 0.0 0.0 0.0 -0.0 -1.14906 -0.381212 0.0 0.0 0.0 … -0.0 -0.0 -0.223315 0.0 0.0 -0.0 0.0105824 -0.0333975 0.0300461 0.0 0.0 0.0 0.11075 0.024764 -0.00163322 0.0 0.0 0.0 0.0 0.0858604 0.028485 0.0 0.0 0.0 0.0 0.0 0.0166866
julia> logdet(C)-17.7305935734267
julia> inv(C)40×40 Kronecker.KroneckerProduct{Float64, Matrix{Float64}, Matrix{Float64}}: 654.355 -307.4 1098.54 … 142.946 -510.842 1403.91 -307.4 536.906 -626.958 -249.671 291.546 -683.213 1098.54 -626.958 2442.4 291.546 -1135.76 2548.79 -3019.06 1469.22 -5481.07 -683.213 2548.79 -7682.66 380.624 -178.808 638.997 83.4225 -298.123 819.313 -178.808 312.306 -364.687 … -145.706 170.144 -398.718 638.997 -364.687 1420.69 170.144 -662.821 1487.46 -1756.12 854.612 -3188.21 -398.718 1487.46 -4483.54 -799.803 375.728 -1342.72 -174.606 623.982 -1714.85 375.728 -656.248 766.316 304.967 -356.117 834.53 ⋮ ⋱ 746.404 -363.237 1355.09 169.098 -630.836 1901.49 -891.688 418.893 -1496.98 -194.875 696.415 -1913.91 418.893 -731.64 854.353 340.369 -397.456 931.404 -1496.98 854.353 -3328.25 -397.456 1548.35 -3474.69 4114.06 -2002.1 7469.03 … 931.404 -3474.69 10473.5 -304.287 142.946 -510.842 -66.8232 238.803 -656.288 142.946 -249.671 291.546 116.714 -136.289 319.382 -510.842 291.546 -1135.76 -136.289 530.934 -1191.48 1403.91 -683.213 2548.79 319.382 -1191.48 3591.41
LinearAlgebra.choleskyFunction
cholesky(K::AbstractKroneckerProduct; check = true)

Wrapper around cholesky from the LinearAlgebra package. Performs Cholesky on the matrices of a AbstractKroneckerProduct instances and returns a CholeskyKronecker type. Similar to Cholesky, size, \, inv, det, and logdet are overloaded to efficiently work with this type.

source