GitHub

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
julia> K = As ⊗ Bs40×40 Kronecker.KroneckerProduct{Float64, Matrix{Float64}, Matrix{Float64}}:
 3.0652      4.42582    0.104919   …   3.42937    0.0812969   4.42638
 4.42582    10.5299    -7.63913        8.15916   -5.91921     0.48486
 0.104919   -7.63913   26.5352        -5.91921   20.5609      9.57074
 5.71254     0.625744  12.3517         0.48486    9.57074    33.6686
 2.56699     3.70646    0.0878656      3.69637    0.0876264   4.77101
 3.70646     8.81841   -6.39748    …   8.7944    -6.38006     0.52261
 0.0878656  -6.39748   22.2222        -6.38006   22.1617     10.3159
 4.78403     0.524037  10.3441         0.52261   10.3159     36.29
 2.40769     3.47645    0.0824129      2.76073    0.0654462   3.56336
 3.47645     8.27117   -6.00047        6.56834   -4.76513     0.390326
 ⋮                                 ⋱
 4.08958     0.447968   8.84251        0.458563   9.05165    31.8426
 2.25335     3.2536     0.0771301      3.16335    0.0749907   4.08303
 3.2536      7.74097   -5.61583        7.52625   -5.46006     0.44725
 0.0771301  -5.61583   19.5071        -5.46006   18.966       8.82834
 4.19951     0.460009   9.0802     …   0.44725    8.82834    31.057
 2.37508     3.42937    0.0812969      4.82567    0.114398    6.22863
 3.42937     8.15916   -5.91921       11.4812    -8.32928     0.682276
 0.0812969  -5.91921   20.5609        -8.32928   28.9325     13.4676
 4.42638     0.48486    9.57074        0.682276  13.4676     47.3772
julia> E = eigen(K)Eigen{Float64, Float64, Kronecker.KroneckerProduct{Float64, Matrix{Float64}, Matrix{Float64}}, Vector{Float64}}
values:
40-element Vector{Float64}:
   0.0035865615341657297
   0.10294632808953044
   0.31850466795597404
   0.6544798540106092
   0.0071723711532584436
   0.20587107369816243
   0.6369425620790219
   1.308822497698036
   0.01819736024426883
   0.5223251853409611
   ⋮
  17.34779154934845
   0.1471583865409296
   4.223938554421704
  13.068403426419462
  26.853630816666744
   1.8420217510569175
  52.872193527404924
 163.58077802352162
 336.13423755087604
vectors:
40×40 Kronecker.KroneckerProduct{Float64, Matrix{Float64}, Matrix{Float64}}:
  0.265047    -0.121324   -0.0612069  …   0.14494     0.0731208   0.0357836
 -0.133933    -0.224156   -0.145105       0.267788    0.17335    -0.0230755
 -0.0282414   -0.139487    0.222866       0.166638   -0.266246    0.167633
 -0.0250488    0.0720452  -0.123053      -0.0860687   0.147005    0.313018
  0.155609    -0.0712299  -0.0359347      0.155845    0.0786219   0.0384758
 -0.0786325   -0.131603   -0.0851915  …   0.287935    0.186391   -0.0248115
 -0.0165806   -0.081893    0.130845       0.179175   -0.286277    0.180244
 -0.0147062    0.0422979  -0.0722444     -0.092544    0.158064    0.336568
 -0.561182     0.25688     0.129593       0.126481    0.0638082   0.0312263
  0.283576     0.474606    0.307231       0.233683    0.151272   -0.0201366
  ⋮                                   ⋱
  0.0325108   -0.0935073   0.15971       -0.0781773   0.133526    0.284319
  0.0908674   -0.0415944  -0.0209839      0.131671    0.0664265   0.0325076
 -0.0459171   -0.0768488  -0.0497472      0.243272    0.157479   -0.0209629
 -0.00968215  -0.0478211   0.0764063      0.151382   -0.241871    0.152286
 -0.00858762   0.0246997  -0.0421868  …  -0.0781889   0.133546    0.284361
 -0.0612208    0.0280237   0.0141376      0.136713    0.0689701   0.0337524
  0.030936     0.0517759   0.0335166      0.252587    0.16351    -0.0217656
  0.00652323   0.0322188  -0.0514777      0.157179   -0.251133    0.158117
  0.00578579  -0.0166411   0.0284228     -0.0811831   0.13866     0.29525
julia> logdet(E)8.34161872675984
julia> b = randn(40);
julia> (E + 0.1I) \ b  # solve a system40-element Vector{Float64}:
 -5.253583408230181
  1.0426467981028615
 -0.5447239178082354
  1.1457778986925373
 -4.156260341817804
  1.2354331914156347
  0.04222846015377739
  0.6559412680247766
  1.2308531028658565
  2.3763084960168666
  ⋮
 -1.0240006580184295
  4.770283130065554
 -2.068765487745894
 -0.4269054173661324
 -0.30087004326909067
  6.76870585331524
 -3.0739676943029455
 -0.5555575169258163
 -0.6487733362643137
LinearAlgebra.eigen — Function
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
Base.:+ — Method
+(E::Eigen, B::UniformScaling)

Shift a spectral decomposition, i.e. eigen(K) + λI yields the (lazy) eigenvalue decomposition of K + λI. This can be used to efficiently solve systems of the form (A ⊗ B + λI) \ v.

source
LinearAlgebra.logdet — Method
logdet(K::Eigen)

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

source
Base.inv — Method
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
julia> K = As ⊗ Bs40×40 Kronecker.KroneckerProduct{Float64, Matrix{Float64}, Matrix{Float64}}:
 17.0938    4.39541  -5.23937    2.77963   …   4.15041  -4.94733    2.62469
  4.39541   6.63697  -5.25932   -4.2224        6.26703  -4.96616   -3.98705
 -5.23937  -5.25932  20.3189    -1.03636      -4.96616  19.1863    -0.978593
  2.77963  -4.2224   -1.03636   14.2686       -3.98705  -0.978593  13.4732
  9.78052   2.51491  -2.99779    1.59041       3.20088  -3.81548    2.02421
  2.51491   3.79745  -3.0092    -2.41592   …   4.83326  -3.83001   -3.07489
 -2.99779  -3.0092   11.6258    -0.59297      -3.83001  14.7969    -0.754711
  1.59041  -2.41592  -0.59297    8.16398      -3.07489  -0.754711  10.3908
  9.8102    2.52254  -3.00689    1.59523       2.02873  -2.41826    1.28295
  2.52254   3.80898  -3.01833   -2.42325       3.06333  -2.42747   -1.94887
  ⋮                                        ⋱
  2.29166  -3.48115  -0.854425  11.7637       -3.71386  -0.911541  12.55
 14.7814    3.8008   -4.53058    2.40359       4.51449  -5.38132    2.85493
  3.8008    5.73911  -4.54783   -3.65119       6.81678  -5.4018    -4.3368
 -4.53058  -4.54783  17.5701    -0.89616      -5.4018   20.8694    -1.06444
  2.40359  -3.65119  -0.89616   12.3383    …  -4.3368   -1.06444   14.6551
 16.141     4.15041  -4.94733    2.62469       5.32737  -6.35027    3.36899
  4.15041   6.26703  -4.96616   -3.98705       8.0442   -6.37444   -5.11768
 -4.94733  -4.96616  19.1863    -0.978593     -6.37444  24.6271    -1.2561
  2.62469  -3.98705  -0.978593  13.4732       -5.11768  -1.2561    17.2939
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}}}:
 4.13447  1.06311  -1.26724   0.672305  2.3656   …  -1.19661     0.634831
 0.0      2.34665  -1.6671   -2.10391   0.0         -1.57417    -1.98664
 0.0      0.0       3.99172  -0.924869  0.0          3.76922    -0.873316
 0.0      0.0       0.0       2.92143   0.0          0.0         2.75858
 0.0      0.0      -0.0       0.0       1.71943     -0.572742    0.303855
 0.0      0.0      -0.0      -0.0       0.0      …  -0.753461   -0.950883
 0.0      0.0       0.0      -0.0       0.0          1.8041     -0.418003
 0.0      0.0       0.0       0.0       0.0          0.0         1.32037
 0.0      0.0      -0.0       0.0       0.0         -0.0750376   0.0398095
 0.0      0.0      -0.0      -0.0       0.0         -0.0987145  -0.12458
 ⋮                                               ⋱
 0.0      0.0       0.0       0.0       0.0         -0.0        -0.709561
 0.0      0.0      -0.0       0.0       0.0         -0.010451    0.00554452
 0.0      0.0      -0.0      -0.0       0.0         -0.0137486  -0.017351
 0.0      0.0       0.0      -0.0       0.0          0.0329198  -0.00762741
 0.0      0.0       0.0       0.0       0.0      …   0.0         0.0240931
 0.0      0.0      -0.0       0.0       0.0         -0.0104021   0.00551858
 0.0      0.0      -0.0      -0.0       0.0         -0.0136843  -0.0172698
 0.0      0.0       0.0      -0.0       0.0          0.0327657  -0.00759173
 0.0      0.0       0.0       0.0       0.0          0.0         0.0239803
julia> logdet(C)-14.296911244817814
julia> inv(C)40×40 Kronecker.KroneckerProduct{Float64, Matrix{Float64}, Matrix{Float64}}:
  1210.72       -1203.13      -29.5227    …    -30.4221     -612.133
 -1203.13        4802.21     1020.98          1052.08       1782.32
   -29.5227      1020.98      960.554          989.817       389.156
  -594.036       1729.62      377.651          389.156      1679.59
 -3582.96        3560.51       87.3687          90.296      1816.87
  3560.51      -14211.5     -3021.46      …  -3122.69      -5290.1
    87.3687     -3021.46    -2842.64         -2937.88      -1155.06
  1757.97       -5118.6     -1117.61         -1155.06      -4985.2
 -1910.57        1898.6        46.5883          48.0211      966.248
  1898.6        -7578.13    -1611.16         -1660.71      -2813.37
     ⋮                                    ⋱
  -735.28        2140.88      467.445          483.898      2088.5
   -29.2921        29.1085      0.714272         0.732411     14.7371
    29.1085      -116.185     -24.7016         -25.3289      -42.9092
     0.714272     -24.7016    -23.2396         -23.8298       -9.36892
    14.3721       -41.8465     -9.13688   …     -9.36892     -40.4361
  1247.6        -1239.78      -30.4221         -31.4975     -633.771
 -1239.78        4948.51     1052.08          1089.27       1845.32
   -30.4221      1052.08      989.817         1024.8         402.912
  -612.133       1782.32      389.156          402.912      1738.96
LinearAlgebra.cholesky — Function
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