Last active
December 29, 2023 11:11
-
-
Save tholden/58dd9a8991daa750ae36a633fe7060a4 to your computer and use it in GitHub Desktop.
Matlab Code for finding the nearest Kronecker product to a matrix
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
function [ B, C, D ] = NearestKroneckerProduct( A, SizeB, SizeC, Hermitian ) %#codegen | |
% https://doi.org/10.1007/978-94-015-8196-7_17 | |
m = size( A, 1 ); | |
n = size( A, 2 ); | |
m1 = SizeB( 1 ); | |
n1 = SizeB( 2 ); | |
m2 = SizeC( 1 ); | |
n2 = SizeC( 2 ); | |
if nargin < 4 | |
Hermitian = false; | |
end | |
assert( m == m1 * m2 ); | |
assert( n == n1 * n2 ); | |
if Hermitian | |
assert( m1 == n1 ); | |
assert( m2 == n2 ); | |
A = 0.5 * ( A + A' ); | |
end | |
R = reshape( permute( reshape( A, [ m2, m1, n2, n1 ] ), [ 2 4 1 3 ] ), m1 * n1, m2 * n2 ); | |
[ B, S, C ] = svds( R, 1 ); | |
SqrtS = sqrt( S ); | |
B = reshape( B * SqrtS, m1, n1 ); | |
C = reshape( C * SqrtS, m2, n2 ); | |
if Hermitian | |
B = 0.5 * ( B + B' ); | |
C = 0.5 * ( C + C' ); | |
if all( diag( B ) < 0 ) && all( diag( C ) < 0 ) | |
B = -B; | |
C = -C; | |
end | |
end | |
if nargout > 2 | |
D = A - kron( B, C ); | |
if Hermitian | |
D = 0.5 * ( D + D' ); | |
end | |
end | |
end |
Testing whether taking Cholesky decompositions first changes things in the exact case:
while true;
m = 2 + poissrnd( 5 );
n = 2 + poissrnd( 5 );
if m==n
continue
end
B = randn( m );
B = B * B';
C = randn( n );
C = C * C';
A = kron( B, C );
[ b, c, d ] = NearestKroneckerProduct( A, [ m m ], [ n n ], true );
[ Lb, Lc, Ld ] = NearestKroneckerProduct( chol( A, 'lower' ), [ m m ], [ n n ], false );
b2 = Lb * Lb';
c2 = Lc * Lc';
if max( max( max( abs( b2 / mean( diag( b2 ) ) - b / mean( diag( b ) ) ) ) ), max( max( abs( c2 / mean( diag( c2 ) ) - c / mean( diag( c ) ) ) ) ) ) > 1e-8
break
end
end
Sign up for free
to join this conversation on GitHub.
Already have an account?
Sign in to comment
Testing whether the decomposition of a p.d. matrix is p.d.: