Edit on 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,Array{Float64,2},Array{Float64,2}}:
 31.7798     13.6217    8.46737   …  11.3318    7.04394    0.0576194
 13.6217     35.7485   13.4123       29.7389   11.1576    -2.68683
  8.46737    13.4123    6.97464      11.1576    5.80216   -1.25406
  0.069263   -3.22978  -1.50748      -2.68683  -1.25406    0.996335
 19.8465      8.50674   5.28787       8.69874   5.40722    0.044231
  8.50674    22.3249    8.37596   …  22.8288    8.56501   -2.06252
  5.28787     8.37596   4.35566       8.56501   4.45397   -0.962666
  0.0432547  -2.017    -0.941418     -2.06252  -0.962666   0.764827
 33.2291     14.2429    8.85351      13.8031    8.58015    0.0701855
 14.2429     37.3788   14.0239       36.2246   13.5909    -3.2728
  ⋮                               ⋱
  0.0503263  -2.34675  -1.09533      -2.24453  -1.04762    0.83232
 23.6713     10.1461    6.30694      10.9539    6.80904    0.0556978
 10.1461     26.6274    9.99016      28.7472   10.7855    -2.59723
  6.30694     9.99016   5.19508      10.7855    5.60866   -1.21224
  0.0515907  -2.40571  -1.12285   …  -2.59723  -1.21224    0.963108
 26.4374     11.3318    7.04394      13.8427    8.60477    0.0703869
 11.3318     29.7389   11.1576       36.3286   13.6299    -3.28219
  7.04394    11.1576    5.80216      13.6299    7.08782   -1.53194
  0.0576194  -2.68683  -1.25406      -3.28219  -1.53194    1.21711
julia> E = eigen(K)Eigen{Float64,Float64,Kronecker.KroneckerProduct{Float64,Array{Float64,2},Array{Float64,2}},Array{Float64,1}}
values:
40-element Array{Float64,1}:
   0.00029123186973506373
   0.0007509239202255859
   0.00991653343638598
   0.025678026641621814
   0.0015982518672834534
   0.00412099664360276
   0.05442096050854867
   0.1409184855540046
   0.03478203104502254
   0.08968338228057506
   ⋮
  22.545207242769884
   0.4543648145901654
   1.1715524406549505
  15.471259654033021
  40.06152151093757
   4.487731491804558
  11.571346665498405
 152.80861752007698
 395.68502208166285
vectors:
40×40 Kronecker.KroneckerProduct{Float64,Array{Float64,2},Array{Float64,2}}:
 -0.00862273    0.00652955  …  -0.0263283  -0.266943    0.20141
 -0.00462438    0.0279036      -0.112512    0.199656    0.246693
  0.0385878    -0.0685141       0.276262    0.0350337   0.109405
  0.0735536     0.0384637      -0.155093   -0.0371207  -0.0182751
  0.00162764   -0.00123253     -0.0209617  -0.212531    0.160355
  0.000872905  -0.00526712  …  -0.0895784   0.158959    0.196408
 -0.00728389    0.0129328       0.21995     0.0278926   0.0871044
 -0.0138841    -0.00726047     -0.123479   -0.0295542  -0.01455
 -0.0508861     0.0385334      -0.0324542  -0.329054    0.248272
 -0.0272903     0.16467        -0.138691    0.24611     0.304091
  ⋮                         ⋱
 -0.357098     -0.186739       -0.129245   -0.0309341  -0.0152293
 -0.0265994     0.0201423      -0.025839   -0.261982    0.197666
 -0.0142653     0.086077       -0.110421    0.195945    0.242108
  0.119036     -0.211352        0.271127    0.0343825   0.107372
  0.226898      0.118653    …  -0.15221    -0.0364308  -0.0179354
 -0.0100118     0.00758142     -0.0261347  -0.26498     0.199928
 -0.00536934    0.0323987      -0.111685    0.198188    0.244878
  0.0448041    -0.0795514       0.27423     0.034776    0.1086
  0.0854027     0.04466        -0.153952   -0.0368477  -0.0181407
julia> logdet(E)-19.21091336330567
julia> b = randn(40);
julia> (E + 0.1I) \ b  # solve a system40-element Array{Float64,1}:
  0.15616820987697325
  1.406985731197416
 -3.991048341428295
 -6.479042238384197
  4.76878233395609
  5.465363034162523
  0.311479955778969
 -9.043735671076762
 -3.2257428410160234
  1.8259458626033038
  ⋮
 -4.4348852139786406
 -2.208750824174579
  2.153671456071821
  1.3938419616001683
  9.190464447632992
  0.7974636152090406
  1.621447536821078
 -2.776927987408376
 -5.341729534106016
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
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.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,Array{Float64,2},Array{Float64,2}}:
  31.9113   -13.1507    5.84077  -2.68873  …   -8.18535   3.63545  -1.67354
 -13.1507    20.5303   -7.48847  -2.72384      12.7786   -4.66102  -1.69539
   5.84077   -7.48847  14.8242    5.42263      -4.66102   9.227     3.37519
  -2.68873   -2.72384   5.42263   8.98161      -1.69539   3.37519   5.59039
  18.4252    -7.59306   3.37239  -1.55244      -7.47186   3.31856  -1.52766
  -7.59306   11.8539   -4.32376  -1.57271  …   11.6647   -4.25474  -1.54761
   3.37239   -4.32376   8.55933   3.13096      -4.25474   8.42271   3.08099
  -1.55244   -1.57271   3.13096   5.18587      -1.54761   3.08099   5.1031
  17.3994    -7.1703    3.18462  -1.46601      -7.29688   3.24084  -1.49189
  -7.1703    11.1939   -4.08302  -1.48515      11.3915   -4.1551   -1.51136
   ⋮                                       ⋱
  -2.36853   -2.39945   4.77684   7.91197      -2.04634   4.07386   6.74762
  16.3562    -6.74041   2.99369  -1.37811      -5.41211   2.40374  -1.10653
  -6.74041   10.5228   -3.83823  -1.39611       8.44912  -3.08184  -1.12098
   2.99369   -3.83823   7.59818   2.77938      -3.08184   6.10084   2.23166
  -1.37811   -1.39611   2.77938   4.60354  …   -1.12098   2.23166   3.69634
  19.8625    -8.18535   3.63545  -1.67354     -10.8129    4.80246  -2.21076
  -8.18535   12.7786   -4.66102  -1.69539      16.8806   -6.15726  -2.23962
   3.63545   -4.66102   9.227     3.37519      -6.15726  12.189     4.45866
  -1.67354   -1.69539   3.37519   5.59039      -2.23962   4.45866   7.38497
julia> C = cholesky(K)40×40 CholeskyKronecker{Cholesky{Float64,Array{Float64,2}},Cholesky{Float64,Array{Float64,2}}}
U factor:
40×40 Kronecker.KroneckerProduct{Float64,UpperTriangular{Float64,Array{Float64,2}},UpperTriangular{Float64,Array{Float64,2}}}:
 5.64901  -2.32797   1.03394  -0.475965  …  -1.44899    0.643555   -0.296253
 0.0       3.88727  -1.30721  -0.985748      2.41954   -0.813645   -0.613556
 0.0       0.0       3.47079   1.33289       0.0        2.16031     0.829624
 0.0       0.0       0.0       2.45087       0.0        0.0         1.52549
 0.0      -0.0       0.0      -0.0          -0.956145   0.424663   -0.195489
 0.0       0.0      -0.0      -0.0       …   1.59658   -0.536901   -0.404868
 0.0       0.0       0.0       0.0           0.0        1.42553     0.547445
 0.0       0.0       0.0       0.0           0.0        0.0         1.00663
 0.0      -0.0       0.0      -0.0          -0.298975   0.132787   -0.0611271
 0.0       0.0      -0.0      -0.0           0.499232  -0.167882   -0.126597
 ⋮                                       ⋱
 0.0       0.0       0.0       0.0           0.0        0.0         0.214224
 0.0      -0.0       0.0      -0.0           0.894184  -0.397144    0.182821
 0.0       0.0      -0.0      -0.0          -1.49312    0.502108    0.378631
 0.0       0.0       0.0       0.0          -0.0       -1.33315    -0.511969
 0.0       0.0       0.0       0.0       …  -0.0       -0.0        -0.941394
 0.0      -0.0       0.0      -0.0          -0.163986   0.0728328  -0.0335278
 0.0       0.0      -0.0      -0.0           0.273825  -0.0920823  -0.0694378
 0.0       0.0       0.0       0.0           0.0        0.244488    0.0938908
 0.0       0.0       0.0       0.0           0.0        0.0         0.172644
julia> logdet(C)18.895749218908648
julia> inv(C)40×40 Kronecker.KroneckerProduct{Float64,Array{Float64,2},Array{Float64,2}}:
  0.843202     0.513028    …   -0.157144    0.0873973   -0.177741
  0.513028     1.3233          -0.405336   -0.103532    -0.107461
 -0.285326     0.338002        -0.103532   -0.557947     0.331625
  0.580271     0.350827        -0.107461    0.331625    -0.863539
  0.00767366   0.00466887      -0.897839    0.499342    -1.01552
  0.00466887   0.0120429   …   -2.31588    -0.591529    -0.613974
 -0.00259664   0.00307603      -0.591529   -3.18782      1.89473
  0.00528082   0.00319274      -0.613974    1.89473     -4.93381
  1.02083      0.6211          -7.76293     4.31744     -8.78042
  0.6211       1.60206        -20.0236     -5.1145      -5.30857
  ⋮                        ⋱
 -0.469413    -0.283803         2.12064    -6.54434     17.0412
 -0.307139    -0.186872         5.53911    -3.08063      6.26512
 -0.186872    -0.482017        14.2875      3.64937      3.78783
  0.103931    -0.123118         3.64937    19.6668     -11.6893
 -0.211366    -0.12779     …    3.78783   -11.6893      30.4385
 -0.258278    -0.157144         6.10541    -3.39559      6.90565
 -0.157144    -0.405336        15.7482      4.02247      4.17509
  0.0873973   -0.103532         4.02247    21.6775     -12.8844
 -0.177741    -0.107461         4.17509   -12.8844      33.5505
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