Source code for spalor.algorithms.mc_algorithms

from math import *

from scipy.sparse import coo_matrix

from spalor.regularization.thresholding import *
from spalor.matrix_tools.factorization_util import *
import numpy as np
from scipy.optimize import minimize

[docs]def lmafit(d1, d2, r, known, data): ''' A rank-constrained matrix completion algorithm that uses a successive over-relaxation scheme. Links for more details: - http://lmafit.blogs.rice.edu/ Parameters ---------- d1 : int number or rows in matrix d2 : int number of columns in matrix r : int target rank of matrix. known : np array array with 2 columns and many rows, indicating indices of the matrix you have a measurement of data : np array vector of measurements, in same order as 'known' References ---------- .. [1] Wen, Z., Yin, W. & Zhang, Y. Solving a low-rank factorization model for matrix completion by a nonlinear successive over-relaxation algorithm. Math. Prog. Comp. 4, 333–361 (2012). https://doi.org/10.1007/s12532-012-0044-1 ''' # set parameters # TODO: parameter selection tol = 1e-5; maxit = 1000; iprint = 2; est_rank = 1; rank_max = max(floor(0.1 * min(d1,d2)), 2 * r); rank_min = 1; rk_jump = 10; init = 0; save_res = 0; # Initialize Variables X = np.random.randn(d1, r) Y = np.random.randn(d2, r) Res = data - partXY(X, Y, known) S = coo_matrix((Res, known), shape=(d1,d2)) alf = 0 increment = 0.1 Res = data - partXY(X, Y, known) res = np.linalg.norm(Res); # main loop for iter in range(0, maxit): X0 = X Y0 = Y Res0 = Res res0 = res X = X + S.dot(Y).dot(np.linalg.pinv(Y.transpose().dot(Y))) XXInv = np.linalg.pinv(X.transpose().dot(X)); Y = Y.dot(X0.transpose().dot(X)).dot(XXInv) + S.transpose().dot(X).dot(XXInv) Res = data - partXY(X, Y, known) res = np.linalg.norm(Res); ratio = res / res0; ''' print("iter: ",iter," residual: ",res,"alf", alf) if ratio >= 1: increment = max(0.1 * alf, 0.1 * increment) X = X0 Y = Y0 Res = Res0 res = res0 alf = 0 elif ratio > 0.7: increment = max(increment, 0.25 * alf) alf = alf + increment ''' S = coo_matrix(((alf + 1) * Res, known)) return (X, Y)
[docs]def alt_proj(m,n,r,X,y): ''' A very simple matrix completion algorithm comprising of two steps at each iteration: - project L onto set of matrices satisfying the measurements - project L onto the set of rank r matrices Parameters ---------- m : int number or rows in matrix n : int number of columns in matrix r : int target rank of matrix. X : np array array with 2 columns and many rows, indicating indices of the matrix you have a measurement of y : np array vector of measurements, in same order as 'known' References ---------- .. [2] Tanner, J., & Wei, K. (2013). Normalized iterative hard thresholding for matrix completion. SIAM Journal on Scientific Computing, 35(5), S104-S125. ''' L=np.zeros((m,n)) for iter in range(0,2000): L[X[0,:],X[1,:]]=y; L=lowRankProj(L,r+max(0,round(10-iter/100))); # print(np.linalg.norm(L[X[0,:],X[1,:]]-y)) u,s,vt=svds(L,r) U=u.dot(np.sqrt(s)); V=np.sqrt(s).dot(vt); return (U,V)
[docs]def svt(m, n, beta_max, known, data, eps=1e-5, r_max=None): ''' Singular value thresholding for matrix completion Parameters ---------- m : int number or rows in matrix n : int number of columns in matrix beta_max : float largest singular value to keep. Larger values mean less regularization, and less estimtor bias known : np array array with 2 columns and many rows, indicating indices of the matrix you have a measurement of data : np array vector of measurements, in same order as 'known' eps : float stopping criteria. (default: 1e-5) r_max : int upper bound on the rank of the matrix. The smaller this is, the faster the algorithm will be References ---------- .. [3] Cai, J. F., Candès, E. J., & Shen, Z. (2010). A singular value thresholding algorithm for matrix completion. SIAM Journal on optimization, 20(4), 1956-1982. ''' if r_max is None: r_max = min(m, n) beta = beta_max / (1.2 ** 30) maxIter = 100; X = np.zeros((m, n)) for iter in range(0, maxIter): X[known] = data X = lowRankSoftThresholding(X, 1 / beta, r_max) beta = min(beta_max, beta * 1.2) return X
[docs]def alt_min(m,n,r, Omega, known): ''' A very simple algorithm for matrix completion via alternating minimization. Parameters ---------- m : int number or rows in matrix n : int number of columns in matrix Omega : np array array with 2 columns and many rows, indicating indices of the matrix you have a measurement of known : np array vector of measurements, in same order as 'known' ''' U=np.random.rand(m,r) V=np.random.rand(r,n) for i in range(0,100): objU=lambda x: np.linalg.norm(np.reshape(x, [m,r]).dot(V)[Omega]-known)**2 U = np.reshape(minimize(objU, U).x, [m,r]) objV=lambda x: np.linalg.norm(U.dot(np.reshape(x, [r,n]))[Omega]-known)**2 V = np.reshape(minimize(objV, V).x, [r,n]) res=np.linalg.norm(U.dot(V)[Omega]-known) if res < 0.0001: break return (U,V)