In [4]:

    

import numpy as np
import matplotlib.pyplot as pl
from tools.plot import matshow
from tools.helpers import Progress, Timer
from physics.ccg_haar import orthogonal
import cvxpy as cvx
import scipy.sparse as sp


    

In [5]:

    

def random_lowrank_hmatrix(n, r):
    A = np.random.uniform(size=(n, r))
    return np.dot(A, A.T)

def random_lowrank_matrix(n, r):
    A, B = np.random.uniform(size=(2, n, r))
    return np.dot(A, B.T)


    

In [6]:

    

def norm_F(A):
    return np.sqrt(np.sum(A**2))


    

In [7]:

    

def choices(n, m):
    res = set()
    while len(res) < m:
        res.add(tuple(np.random.randint(0, n, 2)))
    return res


    

In [8]:

    

def complete_nucnorm(obs, n):
    Z = cvx.Variable(n, n)
    aim = cvx.Minimize(cvx.normNuc(Z))
    cons = [Z[i, j] == val for i, j, val in izip(obs.row, obs.col, obs.data)]
    prob = cvx.Problem(aim, cons)
    try:
        prob.solve(solver=cvx.MOSEK)
    except cvx.SolverError:
        return np.zeros(n, n)
    return Z.value


    

In [9]:

    

def truncate_coo(mat, m):
    return sp.coo_matrix((mat.data[:m], (mat.row[:m], mat.col[:m])),
                        shape=mat.shape)


    

In [10]:

    

import itertools as it
from numpy.linalg import lstsq

def complete_altmin(obs, n, r, regc=1.):
    selmat = np.zeros((n, n), dtype=bool)
    obs_csr = obs.tocsr()
    obs_csc = obs.tocsc()

    L = np.random.randn(n, r)
    R = np.random.randn(n, r)

    while True:
        for i in range(n):
            if obs_csr.indptr[i] == obs_csr.indptr[i + 1]:
                R[i, :] = 0
            else:
                sel = slice(obs_csr.indptr[i], obs_csr.indptr[i + 1])
                L_sel = L[obs_csr.indices[sel]]
                tmat = np.tensordot(L_sel, L_sel, axes=(0, 0)) + regc * np.eye(r)
                sumterm = np.tensordot(L_sel, obs_csr.data[sel], axes=(0, 0))
                R[i, :] = np.linalg.pinv(tmat).dot(sumterm)

        for i in range(n):
            if obs_csc.indptr[i] == obs_csc.indptr[i + 1]:
                L[i, :] = 0
            else:
                sel = slice(obs_csc.indptr[i], obs_csc.indptr[i + 1])
                R_sel = L[obs_csc.indices[sel]]
                tmat = np.tensordot(R_sel, R_sel, axes=(0, 0)) + regc * np.eye(r)
                sumterm = np.tensordot(R_sel, obs_csc.data[sel], axes=(0, 0))
                L[i, :] = np.linalg.pinv(tmat).dot(sumterm)
        
        yield np.dot(L, R.T)


    

In [14]:

    

def iterative_cauchy(iterator, distfunc, reldist=1e-3, maxiter=1000):
    last = next(iterator)
    for _ in range(maxiter):
        current = next(iterator)
        if distfunc(last, current) < reldist:
            return current
        last = current
    return current


    

In [15]:

    

n = 500
m_max = 1000
r = 10

X = random_lowrank_matrix(n, r)
indices = np.array(list(choices(n, m_max)))
data = [X[i, j] for i, j in indices]
obs = sp.coo_matrix((data, (indices[:, 0], indices[:, 1])),
                   shape=(n, n))


    

In [ ]:

    

with Timer():
    solutions = (iterative_cauchy(complete_altmin(truncate_coo(obs, m), n, 2 * r, regc=1. / n),
                                  lambda x,y: norm_F(x - y),
                                  1e-4, maxiter=10000)
                 for m in Progress(range(1, m_max, 100)))
    diffs_am = [norm_F(X - recons) for recons in solutions]


    

    




 20% ( 2 of 10) |##########################                                                                                                          | Elapsed Time: 0:00:00 ETA:  0:00:01

In [ ]:

    

with Timer():
    diffs_no = [norm_F(X - complete_nucnorm(truncate_coo(obs, m), n)) 
                for m in Progress(range(1, m_max, 5))]


    

In [ ]:

    

pl.plot(diffs_am)
#pl.plot(diffs_no)


    

In [ ]:

    

*