##################################### # matrixCompletion_lowrank.mp # # # # By: Pierre-Jean Spaenlehauer # # # # Created: September 19th, 2014 # # Last Revision: # ##################################### ###################################################################### # MatrixCompletion ###################################################################### # Input: - M a matrix with real coefficients and variables (for the unknown entries) # - target_rank # - vars the list of variables # - start: the list of the starting values for the variables # - eps terminating criterion: stop iterating when the size of # the correction step is < eps # - nbitmax: stop if it did not converge after nbitmax # iterations ###################################################################### # Output: - N the completed matrix # - the perturbation ||N-M'|| where M' is the evaluation # of M at the starting point # - number of iterations performed # ###################################################################### MatrixCompletion:=proc(M, target_rank, samples_pos, eps, nbitmax) local Ma,p,q,st,distit,it,timeit,Mold,M2,Bim,Bker,sigma,Tdet,BBim,BBker,i,j,k,A,b,S; Ma:=Matrix(M); p:=LinearAlgebra[RowDimension](M); q:=LinearAlgebra[ColumnDimension](M); st:=time(); distit:=2*eps^2; it:=0; while distit>eps^2 and it<=nbitmax do timeit:=time(); Mold:=Matrix(Ma); userinfo(1,MatrixCompletion,`---------------------------------------`); userinfo(1,MatrixCompletion,`Iteration no:`,it); st:=time(); M2, Bim,Bker,sigma:=truncRank_variant(Ma,target_rank): userinfo(1,MatrixCompletion,`1. Truncating rank. Time:`,time()-st); userinfo(1,MatrixCompletion,`r+1 singular value:`,sigma); st:=time(); Tdet:=[]; BBim:=[seq(LinearAlgebra[SubMatrix](Bim,[1..p],[i..i]),i=1..p)]; BBker:=[seq(LinearAlgebra[Transpose](LinearAlgebra[SubMatrix](Bker,[1..q],[j..j])),j=1..q)]; for i from target_rank+1 to p do for j from 1 to target_rank do Tdet:=[op(Tdet), BBim[i].BBker[j]]; od; od; for i from 1 to target_rank do for j from target_rank+1 to q do Tdet:=[op(Tdet), BBim[i].BBker[j]]; od; od; for i from 1 to target_rank do for j from 1 to target_rank do Tdet:=[op(Tdet), BBim[i].BBker[j]]; od; od; userinfo(1,MatrixCompletion, `2. tangent space of the determinantal variety. Time:`,time()-st); st:=time(); A:=Matrix(nops(samples_pos),(p+q-target_rank)*(target_rank)); b:=Matrix(nops(samples_pos),1); for i from 1 to nops(samples_pos) do b[i,1]:=-M2[samples_pos[i][1],samples_pos[i][2]]+Ma[samples_pos[i][1],samples_pos[i][2]]; for j from 1 to (p+q-target_rank)*target_rank do A[i,j]:=Tdet[j][samples_pos[i][1],samples_pos[i][2]]; od; od; A; b; userinfo(1,MatrixCompletion,`3. Construct linear system. Time:`,time()-st); st:=time(); S:=LinearAlgebra[MatrixInverse](A, method=pseudo).b; userinfo(1,MatrixCompletion,`4. Solving linear system. Time:`,time()-st); st:=time(); it:=it+1; Ma:=M2+add(S[i,1]*Tdet[i],i=1..(p+q-target_rank)*target_rank); for i from 1 to nops(samples_pos) do Ma[samples_pos[i][1],samples_pos[i][2]]:=M[samples_pos[i][1],samples_pos[i][2]]; od; distit:=LinearAlgebra[Trace]((Ma-Mold).LinearAlgebra[Transpose](Ma-Mold)); userinfo(1,MatrixCompletion,`New distance:`,sqrt(distit)); od; Ma,sqrt(LinearAlgebra[Trace]((Ma-M).LinearAlgebra[Transpose](Ma-M))),it; end proc; truncRank_variant:=proc(MM,r) # return the rank r matrix closest to MM and the r+1 th singular value of MM description "rank truncation with SVD (Eckart-Young theorem)"; local U,V,i,S2,ll,m,c, du, dv,S,mm; U,S,V:=MTM[svd](MM): m:=sort([seq(S[i,i], i=1..min([LinearAlgebra[RowDimension](S),LinearAlgebra[ColumnDimension](S)]))]); mm:=m[nops(m)-r+1]; for i from 1 to min(LinearAlgebra[ColumnDimension](S),LinearAlgebra[RowDimension](S)) do if S[i,i] < mm then S[i,i]:=0; fi; od; du:=LinearAlgebra[RowDimension](U); dv:=LinearAlgebra[RowDimension](V); U.S.LinearAlgebra[Transpose](V), LinearAlgebra[SubMatrix](U,[1..du],[1..du]),LinearAlgebra[SubMatrix](V,[1..dv],[1..dv]),m[nops(m)-r]; end proc; # with(LinearAlgebra): # M:=Matrix(3,3,[0,15.4,90.7,-31.3,0,-75.1,56.8,2.4,0]); # spos:=[[1,2],[1,3],[2,1],[2,3],[3,1],[3,2]]; # infolevel[MatrixCompletion]:=1; # MatrixCompletion_variant(M,2,spos,0.0001,50);