PSD Factorization

This code allows you to compute (approximate) PSD factorizations:

Given an m-by-n nonnegative matrix X and an integer k, it computes symmetric k-by-k positive semidefinite matrices {A¹,...,A^m} and {B¹,...,Bⁿ} such that the quantity ∑_i,j ( X(i,j) - trace( Aⁱ B^j ) )² is minimized.

Note that if Aⁱ and B^j are diagonal matrices for all i,j, the above problem reduces to NMF.

See A. Vandaele, F. Glineur and N. Gillis, "Algorithms for Positive Semidefinite Factorization", Computational Optimization and Applications 71 (1), pp. 193-219, 2018. [doi] [arXiv]

Numerical example on the slack matrix of the square:

>> addpath matrices/;

>> X=slackregularngon(4);

>> X

X =

0 1.0000 1.0000 0

0 0 1.0000 1.0000

1.0000 0 0 1.0000

1.0000 1.0000 0 0

% You can compute an (almost) exact PSD factorization of rank 3 (since rank_psd(X)=3):

>> [m,n]=size(X);

>> k = 3;

>> maxiter = Inf;

>> timelimit = 2;

>> [a0,b0,A0,B0] = initfactors(m,n,k);

>> [a,b,e,t] = PSDFact_CD(X,a0,b0,maxiter,timelimit,'greedy');

>> fprintf('Relative error found with CD_gausssouthwell is %1.3g%%\n',100*e(end)/norm(X,'fro'));

Relative error found with CD_gausssouthwell is 0.000226%

Using plot( t, log10(e) ) gives

% You cannot compute an exact PSD factorization of rank 2:

>> k = 2;

>> [a0,b0,A0,B0] = initfactors(m,n,k);

>> [a,b,e,t] = PSDFact_CD(X,a0,b0,maxiter,timelimit,'greedy');

>> fprintf('Relative error found with CD_gausssouthwell is %1.3g%%\n',100*e(end)/norm(X,'fro'));

Relative error found with CD_gausssouthwell is 16.9%