Skip to content

Instantly share code, notes, and snippets.

Embed
What would you like to do?
Matlab Code for finding the nearest Kronecker product to a matrix
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
@tholden

This comment has been minimized.

Copy link
Owner Author

commented Jun 22, 2018

Testing whether the decomposition of a p.d. matrix is p.d.:

while true;
    m = 2 + poissrnd( 5 );
    n = 2 + poissrnd( 5 );
    A = randn( m * n );
    A = A * A';
    [ b, c, d ] = NearestKroneckerProduct( A, [ m m ], [ n n ], true );
    if min( [ eig( b ); eig( c ) ] ) < 0
        break
    end
end
@tholden

This comment has been minimized.

Copy link
Owner Author

commented Jun 22, 2018

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
You can’t perform that action at this time.