Instantly share code, notes, and snippets.

# tholden/NearestKroneckerProduct.m Last active Jun 25, 2018

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

### tholden commented Jun 22, 2018 • edited

 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 ``````
Owner Author

### tholden 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 ``````
to join this conversation on GitHub. Already have an account? Sign in to comment
You can’t perform that action at this time.