Matlab Code for finding the nearest Kronecker product to a matrix

NearestKroneckerProduct.m

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 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

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