% % ---------------------------------------------------------------------------- % Given a matrix A (of size m1*m2, n1*n2), find the optimal kronecker % decomposition of the matrix. % the following procedure finds the B,C pair that minimize % norm( A - kron(B,C), 'fro' )**2 % % See Van Loan and Pitsianis, "Approximation with Kronecker Products." % function [B,C] = krondecomp( A, m1, m2, n1, n2 ) tmpA = rearrange( A, m1, m2, n1, n2 ); [u,a,v] = svd( tmpA ); vecB = a(1,1) * u(:,1); vecC = v(:,1); B = vtm( vecB, m1, n1 ); C = vtm( vecC, m2, n2 ); %norm( A - kron(B,C), 'fro' ) % % ---------------------------------------------------------------------------- % The rearrange operator rearranges a matrix in a form that is suitable % for a rank-one approximation through kronecker products. % % ---------------------------------------------------------------------------- function [newA] = rearrange( A, m1, m2, n1, n2 ) % there should be m1*m2 rows and n1*n2 columns [ h,w ] = size( A ); if ( h ~= m1*m2 ) error( 'Incorrect blocking parameters m1 m2!' ); end; if ( w ~= n1*n2 ) error( 'Incorrect blocking parameters n1 n2!' ); end; for j = 1:n1 Aj = []; for i = 1:m1 tmp = mblock( A, i, j, m2, n2 ); tmpv = mtv( tmp )'; Aj( i, 1:length(tmpv) ) = tmpv; end; [ th, tw ] = size( Aj ); newA( (j-1)*th+1:j*th , 1:tw ) = Aj; end; % ---------------------------------------------------------------------------- % get the i,j block of A (which is of size m2xn2) function [rval] = mblock( A, i, j, m2, n2 ) rval = A( (i-1)*m2 + 1:i*m2, (j-1)*n2+1:j*n2 ); % ---------------------------------------------------------------------------- % converts a matrix to a vector by stacking the columns function [vector] = mtv( A ) [ h,w ] = size( A ); for ii = 1:w vector( (ii-1)*h+1:ii*h ) = A( : , ii ); end vector = vector'; % ---------------------------------------------------------------------------- % converts a vector to a matrix by unstacking the columns function [newA] = vtm( vector, h, w ) l = length( vector ); for i = 1:l/h newA( 1:h, i ) = vector( (i-1)*h+1:i*h ); end % ----------------------------------------------------------------------------