Alternating Least Square


In [2]:
import numpy as np
import pandas as pd
from matplotlib import pylab as plt
import time
%matplotlib inline

In [3]:
## Y = A*B

M = 3
N = 5
K = 2

A_true = np.mat(np.random.randn(M, K))
B_true = np.mat(np.random.randn(K, N))
print(A_true)
print(B_true)


[[-1.21484324 -0.04838232]
 [-0.28746379 -1.26890843]
 [-2.33077118 -0.86539365]]
[[ 2.5809738  -0.36618532 -1.59245333  0.52606297 -0.45481611]
 [-0.73261514 -0.03034256 -0.94484371 -0.32557128  0.31137637]]

In [4]:
Y = A_true * B_true
Y


Out[4]:
matrix([[-3.10003294,  0.4463258 ,  1.98029488, -0.62333215,  0.53746517],
        [ 0.18768502,  0.14376695,  1.65669281,  0.26189609, -0.26436493],
        [-5.38165887,  0.87975244,  4.52930607, -0.9443851 ,  0.79060916]])

In [5]:
B = A_true.I * Y
B


Out[5]:
matrix([[ 2.5809738 , -0.36618532, -1.59245333,  0.52606297, -0.45481611],
        [-0.73261514, -0.03034256, -0.94484371, -0.32557128,  0.31137637]])

In [6]:


In [7]:


In [8]:
U * M


Out[8]:
matrix([[ 2.81930616, -1.16777074, -0.75075461,  4.11140669],
        [-2.23243292,  1.67940427, -1.45611726,  2.4167917 ],
        [-0.47924649,  0.7240146 , -0.42121213,  2.9992647 ],
        [-0.26453334, -0.03093931, -0.10407464, -2.18723647]])

In [9]:
U


Out[9]:
matrix([[ 0.61990172,  1.14702424,  0.72399525, -0.83033268,  0.49684847,
         -1.19278353,  0.24295516,  0.08382553,  0.02033281, -0.93687543],
        [ 0.5324797 , -0.3902074 , -0.79594402, -1.20434216, -0.75490146,
          0.55632023,  1.66170454, -1.76342203,  0.21643463, -0.24178323],
        [-0.99396494,  0.62709088, -0.05774733, -0.78453276, -1.06070693,
         -1.07044867,  0.62747602,  0.53446778, -0.44097795, -0.23960225],
        [ 0.94572044, -0.44020372, -0.73421287, -0.19860962,  1.24212543,
          1.19974624, -0.06604063, -0.81539128, -0.22041883, -0.30599612]])

In [10]:
M


Out[10]:
matrix([[-0.13002397, -0.34584708,  0.48927543,  0.34997509],
        [-0.74952918, -1.1954488 , -0.34901751,  0.98857112],
        [ 1.55863349,  0.53371217, -0.64006659,  0.2568413 ],
        [ 1.19896162, -1.43919556,  0.02838537, -0.089313  ],
        [ 1.48640297, -0.45115731, -0.24884043, -0.62280481],
        [-0.45781105,  0.46390734,  0.19706631, -1.60131805],
        [ 1.71681563, -0.09823042, -0.33815042,  0.9298259 ],
        [ 1.00524998,  0.10167474,  1.2740006 , -0.77658219],
        [ 0.23882929, -0.34020008,  0.66169194,  0.134848  ],
        [-1.95961244,  0.38834324, -0.16421577, -0.78590332]])

In [11]:
U = A * M.I

In [12]:
M = U.I * A

In [13]:
U


Out[13]:
matrix([[ 1.07011144, -0.4519106 , -0.08838942, -1.78034151, -1.99311797,
         -1.74226411,  0.83524173,  1.21926452,  1.04893798, -1.18776432],
        [ 2.77651471, -0.90313191, -0.55161375, -4.07060103, -4.85962699,
         -4.2271474 ,  1.8674498 ,  3.18092834,  2.74723607, -2.88458133],
        [ 2.57404251, -1.35642107,  0.10966667, -4.83110653, -5.10596288,
         -4.48417314,  2.30875886,  2.91539424,  2.49745384, -3.05424027]])

In [14]:
M


Out[14]:
matrix([[-0.14210157, -0.28420314,  0.53678994,  0.25258681],
        [-0.26210432, -0.52420863,  0.29818929, -0.22601934],
        [ 0.39708113,  0.79416226, -0.67224602,  0.12191624],
        [-0.41963405, -0.8392681 ,  0.17595793, -0.66331017],
        [-0.10805474, -0.21610949, -0.39244521, -0.60855469],
        [-0.11938126, -0.23876253, -0.30243527, -0.5411978 ],
        [ 0.24739988,  0.49479976, -0.16488726,  0.32991249],
        [-0.18273735, -0.3654747 ,  0.64554941,  0.28007471],
        [-0.16996895, -0.3399379 ,  0.57615962,  0.23622172],
        [-0.07805965, -0.15611929, -0.21160209, -0.36772138]])

In [15]:
U * M


Out[15]:
matrix([[  1.,   2.,   3.,   5.],
        [  2.,   4.,   8.,  12.],
        [  3.,   6.,   7.,  13.]])

Example: Documents vs Words


In [234]:
A = np.mat([
    # colums: words, rows: documents
 # the, a, sparse, love, war, gun
    [21, 22,  0, 70, 10, 6],
    [22, 21,  0, 70, 10, 6],
    [20, 23,  0,  0, 10, 6],
    [22, 24, 65,  0, 10, 6],
    [21, 22, 65,  0, 10, 6],
])

In [17]:
np.linalg.svd(A)[1]


Out[17]:
array([ 118.04858064,   95.16309661,   26.47771342,    1.10092841,
          0.4861085 ])

In [18]:
np.linalg.svd(A)[1]


Out[18]:
array([ 118.04858064,   95.16309661,   26.47771342,    1.10092841,
          0.4861085 ])

In [19]:
np.linalg.svd(A)[1]


Out[19]:
array([ 118.04858064,   95.16309661,   26.47771342,    1.10092841,
          0.4861085 ])

In [20]:
np.linalg.svd(A)[0][0] np.linalg.svd(A)[1][0]


  File "<ipython-input-20-e7f16e8da997>", line 1
    np.linalg.svd(A)[0][0] np.linalg.svd(A)[1][0]
                            ^
SyntaxError: invalid syntax

In [21]:
np.linalg.svd(A)[0][:,0] * np.linalg.svd(A)[1][0] * np.linalg.svd(A)[2][0]


Out[21]:
matrix([[ 24.98898798,  25.96960455,  30.65397044,  42.46979635,
          11.68439831,   7.01063899],
        [ 24.98319204,  25.96358117,  30.64686056,  42.45994591,
          11.68168824,   7.00901294],
        [  7.42352669,   7.71484034,   9.10643392,  12.61658404,
           3.47110667,   2.082664  ],
        [ 19.64914502,  20.42021575,  24.10358959,  33.39451714,
           9.18758443,   5.51255066],
        [ 19.1944605 ,  19.94768852,  23.54582847,  32.62176239,
           8.97498218,   5.38498931]])

In [22]:
def reconstruct(A, k=1):
    return np.linalg.svd(A)[0][:,:k]  * np.diag(np.linalg.svd(A)[1][:k]) * np.linalg.svd(A)[2][:k]

In [23]:
reconstruct(A, 1)


Out[23]:
matrix([[ 24.98898798,  25.96960455,  30.65397044,  42.46979635,
          11.68439831,   7.01063899],
        [ 24.98319204,  25.96358117,  30.64686056,  42.45994591,
          11.68168824,   7.00901294],
        [  7.42352669,   7.71484034,   9.10643392,  12.61658404,
           3.47110667,   2.082664  ],
        [ 19.64914502,  20.42021575,  24.10358959,  33.39451714,
           9.18758443,   5.51255066],
        [ 19.1944605 ,  19.94768852,  23.54582847,  32.62176239,
           8.97498218,   5.38498931]])

In [24]:
reconstruct(A, 2)


Out[24]:
matrix([[ 22.45010659,  22.65481625,  -0.88700814,  69.20222934,
          10.49270475,   6.29562285],
        [ 22.44317119,  22.64730516,  -0.90827382,  69.20437661,
          10.48945983,   6.2936759 ],
        [  7.63386117,   7.98945509,  11.71945688,  10.40192665,
           3.56983292,   2.14189975],
        [ 22.88259678,  24.64184191,  64.27334036,  -0.65119796,
          10.70529364,   6.42317618],
        [ 22.41444781,  24.15173536,  63.54830798,  -1.28218261,
          10.48637148,   6.29182289]])

In [25]:
reconstruct(A, 3)


Out[25]:
matrix([[  2.15279357e+01,   2.15331217e+01,  -1.24088240e-02,
           6.99784719e+01,   1.00142686e+01,   6.00856115e+00],
        [  2.14724193e+01,   2.14665184e+01,   1.24004347e-02,
           7.00215126e+01,   9.98581907e+00,   5.99149144e+00],
        [  1.99905206e+01,   2.30196393e+01,   2.33961920e-04,
           6.38133216e-04,   9.98065437e+00,   5.98839262e+00],
        [  2.21195642e+01,   2.37137172e+01,   6.49970108e+01,
          -8.91092457e-03,   1.03094208e+01,   6.18565248e+00],
        [  2.08806483e+01,   2.22860783e+01,   6.50029842e+01,
           8.90188988e-03,   9.69061313e+00,   5.81436788e+00]])

SINGULAR VALUE DECOMPOSITION


In [235]:
A = np.mat([
    [5, 0, 5, 3, 0, 0, 1, 0, 0, 1, 0],
    [4, 5, 0, 4, 3, 0, 2, 1, 0, 0, 2],
    [5, 4, 4, 4, 0, 0, 1, 0, 1, 0, 1],
    [0, 1, 0, 0, 0, 3, 0, 0, 0, 3, 0],
    [1, 0, 0, 1, 1, 5, 4, 0, 3, 5, 0],
    [0, 1, 1, 0, 1, 4, 5, 5, 4, 0, 4]
])
shape = A.shape
print("rows, cols: ", shape)


rows, cols:  (6, 11)

In [27]:
U, S, V = np.linalg.svd(A)
(U, S, V)


Out[27]:
(matrix([[-0.37770698,  0.3668459 ,  0.35582859,  0.54372574, -0.25276058,
          -0.48643802],
         [-0.47808234,  0.26586112, -0.34578263, -0.65784554, -0.26685458,
          -0.27788018],
         [-0.49655631,  0.43355096,  0.02103743,  0.091733  ,  0.44317579,
           0.60016922],
         [-0.12090558, -0.18971996,  0.35361987, -0.27074126,  0.727388  ,
          -0.47111255],
         [-0.38879129, -0.49257202,  0.58400825, -0.20947981, -0.35540889,
           0.30814159],
         [-0.46522008, -0.57284785, -0.53597124,  0.38210469,  0.11439913,
          -0.09518034]]),
 array([ 14.5279703 ,  10.68614277,   6.87089724,   5.35101615,
          2.43855306,   1.71909515]),
 matrix([[ -4.59282127e-01,  -3.41600546e-01,  -2.98732730e-01,
           -3.73105583e-01,  -1.57507093e-01,  -2.86864125e-01,
           -3.93151513e-01,  -1.93019580e-01,  -2.42553531e-01,
           -1.84773241e-01,  -2.28084256e-01],
         [  4.27923984e-01,   2.15320130e-01,   2.80324301e-01,
            3.18694415e-01,  -2.50639080e-02,  -4.98160232e-01,
           -3.27752330e-01,  -2.43153977e-01,  -3.12138491e-01,
           -2.49404685e-01,  -1.24096993e-01],
         [  1.57942671e-01,  -2.65920841e-01,   1.93180221e-01,
            5.13052710e-02,  -1.43985689e-01,   2.67361867e-01,
           -9.58422734e-02,  -4.40355707e-01,  -5.39700634e-02,
            6.31173674e-01,  -4.09613573e-01],
         [  6.28724975e-02,  -5.25308124e-01,   6.48038675e-01,
           -1.57493963e-01,  -3.36555093e-01,  -6.18955419e-02,
            7.33266137e-02,   2.34100944e-01,   1.85331594e-01,
           -2.45915363e-01,   5.68977350e-02],
         [ -1.93045285e-01,   5.24990580e-01,   2.55602143e-01,
           -1.67478741e-01,  -4.27127676e-01,   3.53781963e-01,
           -4.89197409e-01,   1.25132025e-01,  -6.78493895e-02,
            6.24792526e-02,   1.50524984e-01],
         [ -1.36538760e-01,   2.58847269e-01,  -7.36978111e-02,
            8.02653042e-02,  -3.61049991e-01,  -1.47374665e-01,
            1.83023905e-01,  -4.38475947e-01,   6.65392275e-01,
           -2.08870192e-01,  -1.95633438e-01],
         [ -3.89792049e-01,   3.07148703e-01,   4.86349420e-01,
           -3.20549454e-01,   3.04787117e-01,  -1.67862457e-01,
            4.13381949e-01,  -2.13775926e-01,  -2.45352969e-01,
            6.54795560e-02,  -1.10863413e-01],
         [  1.62651374e-01,   1.93882703e-01,  -1.01411404e-01,
           -2.63279539e-01,  -1.81299739e-01,  -5.00251057e-01,
            4.80152762e-02,   5.12382606e-01,   1.37668212e-01,
            4.35623490e-01,  -3.15707398e-01],
         [ -1.59317170e-02,  -2.57038701e-02,   1.95997059e-01,
           -1.44610576e-01,   6.22480417e-01,  -4.51250690e-02,
           -5.20188010e-01,   4.04309738e-05,   5.25324711e-01,
            5.36930257e-02,  -2.82085697e-02],
         [ -5.71133597e-01,  -1.03219167e-01,   1.29896404e-01,
            6.65564933e-01,  -4.06604917e-02,  -2.57336895e-01,
           -8.22521187e-02,   1.55757149e-01,   3.44018952e-02,
            2.91743283e-01,   1.34549525e-01],
         [  1.27640600e-01,  -6.64910420e-02,  -4.50837235e-02,
           -2.47670950e-01,  -1.13432941e-01,  -3.17275097e-01,
           -9.21030889e-03,  -3.38360210e-01,   5.31620413e-02,
            3.39438778e-01,   7.54828132e-01]]))

In [28]:
def reconstruct(U, S, V, rank):
    return U[:,:rank]  * np.diag(S[:rank]) * V[:rank]

In [29]:
reconstruct(U, S, V, 1)


Out[29]:
matrix([[ 2.52022607,  1.87447007,  1.63924083,  2.04734815,  0.86429116,
          1.57411404,  2.15734651,  1.05915939,  1.33096782,  1.01390912,
          1.25157034],
        [ 3.18997431,  2.37260913,  2.07486788,  2.59142944,  1.09397591,
          1.9924337 ,  2.73065977,  1.34063023,  1.68467155,  1.28335473,
          1.58417425],
        [ 3.31324076,  2.46429109,  2.15504458,  2.69156702,  1.13624914,
          2.06942499,  2.83617747,  1.39243464,  1.74977035,  1.33294591,
          1.64538964],
        [ 0.80673486,  0.60002568,  0.52472782,  0.65536467,  0.2766632 ,
          0.50388046,  0.6905756 ,  0.33904133,  0.42604834,  0.32455653,
          0.40063288],
        [ 2.59418541,  1.92947886,  1.6873465 ,  2.10743027,  0.88965492,
          1.62030849,  2.22065667,  1.09024181,  1.37002682,  1.04366361,
          1.28829932],
        [ 3.10415171,  2.3087768 ,  2.01904595,  2.52171   ,  1.06454374,
          1.93882956,  2.6571945 ,  1.30456211,  1.63934739,  1.24882754,
          1.54155386]])

In [30]:
def calcError(A, rA):
    return np.sum(np.power(A - rA, 2))

In [31]:
S_sum = np
for i in range(1, len(S)):
    rA = reconstruct(U, S, V, i)
    percentage = S[:i].sum() / S.sum()
    error =  calcError(A, rA)
    print("with rank {}, coverrage: {:0.4f}, error= {:0.4f}".format(i, percentage, error))


with rank 1, coverrage: 0.3493, error= 198.9381
with rank 2, coverrage: 0.6062, error= 84.7444
with rank 3, coverrage: 0.7714, error= 37.5352
with rank 4, coverrage: 0.9000, error= 8.9018
with rank 5, coverrage: 0.9587, error= 2.9553

In [32]:
reconstruct(U, S, V, 4)


Out[32]:
matrix([[ 4.76683441,  0.54004518,  5.09591696,  2.96389163, -0.56519073,
          0.09482103,  0.85152402, -0.28954052,  0.5146028 ,  0.86384615,
         -0.07081608],
        [ 3.80915299,  5.46528385,  0.13112467,  3.92935798,  2.54957687,
          0.1598185 ,  1.76909111,  0.87196726,  0.27370728, -0.05912012,
          2.00449791],
        [ 5.34949911,  3.1655735 ,  3.79980637,  4.098182  ,  0.83411276,
         -0.23028137,  1.33984485,  0.31716553,  0.38680822,  0.14797959,
          1.0391709 ],
        [ 0.23183774,  0.27842236, -0.51306741,  0.36207533,  0.46521826,
          2.2531138 ,  1.01595449, -0.57707206,  0.65924225,  2.7200146 ,
         -0.42543837],
        [ 0.90501881,  0.31788284,  0.26056569,  0.81233044,  0.82107262,
          5.38468497,  3.47906843,  0.3407214 ,  2.58872133,  5.16479343,
          0.23408938],
        [ 0.0315125 ,  0.89589788,  0.91663636,  0.0598546 ,  1.06007864,
          3.87719199,  5.16641771,  4.89334691,  4.12780202, -0.05160592,
          3.92599796]])

In [33]:
A


Out[33]:
matrix([[5, 0, 5, 3, 0, 0, 1, 0, 0, 1, 0],
        [4, 5, 0, 4, 3, 0, 2, 1, 0, 0, 2],
        [5, 4, 4, 4, 0, 0, 1, 0, 1, 0, 1],
        [0, 1, 0, 0, 0, 3, 0, 0, 0, 3, 0],
        [1, 0, 0, 1, 1, 5, 4, 0, 3, 5, 0],
        [0, 1, 1, 0, 1, 4, 5, 5, 4, 0, 4]])

In [34]:
A.shape


Out[34]:
(6, 11)

In [35]:
mask = np.ones(A.shape)
mask[A==0] = 0
mask


Out[35]:
array([[ 1.,  0.,  1.,  1.,  0.,  0.,  1.,  0.,  0.,  1.,  0.],
       [ 1.,  1.,  0.,  1.,  1.,  0.,  1.,  1.,  0.,  0.,  1.],
       [ 1.,  1.,  1.,  1.,  0.,  0.,  1.,  0.,  1.,  0.,  1.],
       [ 0.,  1.,  0.,  0.,  0.,  1.,  0.,  0.,  0.,  1.,  0.],
       [ 1.,  0.,  0.,  1.,  1.,  1.,  1.,  0.,  1.,  1.,  0.],
       [ 0.,  1.,  1.,  0.,  1.,  1.,  1.,  1.,  1.,  0.,  1.]])

In [36]:
A_true = np.mat(np.random.randn(6, 5))
B_true = np.mat(np.random.randn(5, 11))

In [37]:
B_true = A.I * A
A_true = A * B_true.I
calcError(A, A_true * B_true)


Out[37]:
9952.3552777548921

In [38]:
B_true = A.I * A
A_true = A * B_true.I
calcError(A, A_true * B_true)


Out[38]:
9952.3552777548921

In [39]:
A


Out[39]:
matrix([[5, 0, 5, 3, 0, 0, 1, 0, 0, 1, 0],
        [4, 5, 0, 4, 3, 0, 2, 1, 0, 0, 2],
        [5, 4, 4, 4, 0, 0, 1, 0, 1, 0, 1],
        [0, 1, 0, 0, 0, 3, 0, 0, 0, 3, 0],
        [1, 0, 0, 1, 1, 5, 4, 0, 3, 5, 0],
        [0, 1, 1, 0, 1, 4, 5, 5, 4, 0, 4]])

In [224]:


In [269]:
A = np.random.random((10, 20))

In [270]:
rank = 4
U1 = np.mat(np.random.randn(A.shape[0], rank))
M1 = np.mat(np.random.randn(rank, A.shape[1]))
U1.shape, M1.shape


Out[270]:
((10, 4), (4, 20))

In [424]:
U = np.copy(U1)
M = np.copy(M1)
errors = []
start = time.time()
for i in range(200):
    U = np.linalg.solve(M.dot(M.T), M.dot(A.T)).T
    M = np.linalg.solve(U.T.dot(U), U.T.dot(A))
    if i % 50 == 0:
        e = np.sqrt(np.power(A - U.dot(M), 2).sum()) / (A.shape[0] + A.shape[1])
        errors.append(e)
elapsed = time.time() - start
e = np.sqrt(np.power(A - U.dot(M), 2).sum()) / (A.shape[0] + A.shape[1])
errors.append(e)
print(elapsed, errors[-1])
plt.plot(errors)
plt.show()


0.014291763305664062 0.0787419623936

Stochastic Gradient Descent


In [451]:
U = np.copy(U1)
M = np.copy(M1)
eta = 0.2
MAX_ITERATION = 20000
errors = []
start = time.time()
for i in range(MAX_ITERATION):
    if i % 200 == 0:
        eta = eta * 0.9
    for j in range(A.shape[0]):
#         u_ind = np.random.randint(0, A.shape[0])
        u_ind = j
        m_ind = np.random.randint(0, A.shape[1])
        a_i = A[u_ind, m_ind]
        u_i = U[u_ind,:]
        m_i = M[:,m_ind]

        err = a_i - u_i.dot(m_i)    
        dU = - err * m_i.T
        dM = - u_i.T.dot(err)

        U[u_ind,:] = U[u_ind,:] - eta * dU
        M[:,m_ind] = M[:,m_ind] - eta * dM
    if i % 400 == 0:
        e = np.sqrt(np.power(A - U.dot(M), 2).sum()) / (A.shape[0] + A.shape[1])
        errors.append(e)
elapsed = time.time() - start
e = np.sqrt(np.power(A - U.dot(M), 2).sum()) / (A.shape[0] + A.shape[1])
errors.append(e)
print(elapsed, errors[-1])
plt.plot(errors)
plt.show()


4.065443992614746 0.0793941976917

In [452]:
U = np.copy(U1)
M = np.copy(M1)
eta = 0.14
MAX_ITERATION = 40000
block_size = 4
errors = []
start = time.time()
for i in range(MAX_ITERATION):
#     for j in range(0, A.shape[0] - block_size):
#         u_ind = j
    u_ind = np.random.randint(0, A.shape[0]-block_size)
    m_ind = np.random.randint(0, A.shape[1]-block_size)
    a_i = A[u_ind:u_ind+block_size, m_ind:m_ind+block_size]
    u_i = U[u_ind:u_ind+block_size,:]
    m_i = M[:,m_ind:m_ind+block_size]

    err = a_i - u_i.dot(m_i)
    dU = - err.dot(m_i.T)
    dM = - u_i.T.dot(err)

    U[u_ind:u_ind+block_size,:] = U[u_ind:u_ind+block_size,:] - eta * dU
    M[:,m_ind:m_ind+block_size] = M[:,m_ind:m_ind+block_size] - eta * dM
    if i % 100 == 0:
        e = np.sqrt(np.power(A - U.dot(M), 2).sum()) / (A.shape[0] + A.shape[1])
        errors.append(e)
elapsed = time.time() - start
print('k')
e = np.sqrt(np.power(A - U.dot(M), 2).sum()) / (A.shape[0] + A.shape[1])
errors.append(e)
print(elapsed, errors[-1])
plt.plot(errors)
plt.show()


k
1.0532281398773193 0.131717052629

In [458]:
U = np.copy(U1)
M = np.copy(M1)
eta = 0.01
MAX_ITERATION = 12000
block_size = 1
errors = []
start = time.time()
for i in range(MAX_ITERATION):
    for j in range(0, A.shape[0] - block_size):
        u_ind = j
        m_ind = np.random.randint(0, A.shape[1]-block_size)
        a_i = A[u_ind:u_ind+block_size, m_ind:m_ind+block_size]
        u_i = U[u_ind:u_ind+block_size,:]
        m_i = M[:,m_ind:m_ind+block_size]
        
#         u_i = np.linalg.solve(m_i.dot(m_i.T), m_i.dot(a_i.T)).T
        err = a_i - u_i.dot(m_i)
        dU = - err.dot(m_i.T)
        U[u_ind:u_ind+block_size,:] = U[u_ind:u_ind+block_size,:] - eta * dU
#         U[u_ind:u_ind+block_size,:] = u_i
        
    for j in range(0, A.shape[1] - block_size):
        u_ind = np.random.randint(0, A.shape[0]-block_size)
        m_ind = j
        a_i = A[u_ind:u_ind+block_size, m_ind:m_ind+block_size]
        u_i = U[u_ind:u_ind+block_size,:]
        m_i = M[:,m_ind:m_ind+block_size]
#         m_i = np.linalg.solve(u_i.T.dot(u_i), u_i.T.dot(a_i))
        err = a_i - u_i.dot(m_i)
        dM = - u_i.T.dot(err)
        M[:,m_ind:m_ind+block_size] = M[:,m_ind:m_ind+block_size] - eta * dM
        
        
    if i % 2400 == 0:
        e = np.sqrt(np.power(A - U.dot(M), 2).sum()) / (A.shape[0] + A.shape[1])
        errors.append(e)
elapsed = time.time() - start
print('k')
e = np.sqrt(np.power(A - U.dot(M), 2).sum()) / (A.shape[0] + A.shape[1])
errors.append(e)
print(elapsed, errors[-1])
plt.plot(errors)
plt.show()


k
4.646364212036133 0.379032177699

In [462]:
U = np.copy(U1)
M = np.copy(M1)
eta = 0.01
MAX_ITERATION = 12000
block_size = 1
errors = []
start = time.time()
for i in range(MAX_ITERATION):
    for j in range(0, A.shape[0] - block_size):
        u_ind = j
        m_ind = np.random.randint(0, A.shape[1]-block_size)
        a_i = A[u_ind:u_ind+block_size, :]
        u_i = U[u_ind:u_ind+block_size,:]
        
        err = a_i - u_i.dot(M)
        dU = - err.dot(M.T)
        U[u_ind:u_ind+block_size,:] = U[u_ind:u_ind+block_size,:] - eta * dU
#         U[u_ind:u_ind+block_size,:] = u_i
        
    for j in range(0, A.shape[1] - block_size):
        u_ind = np.random.randint(0, A.shape[0]-block_size)
        m_ind = j
        a_i = A[:, m_ind:m_ind+block_size]
        m_i = M[:,m_ind:m_ind+block_size]

        err = a_i - U.dot(U)
        dM = - U.T.dot(err)
        M[:,m_ind:m_ind+block_size] = M[:,m_ind:m_ind+block_size] - eta * dM
        
        
    if i % 2400 == 0:
        e = np.sqrt(np.power(A - U.dot(M), 2).sum()) / (A.shape[0] + A.shape[1])
        errors.append(e)
elapsed = time.time() - start
print('k')
e = np.sqrt(np.power(A - U.dot(M), 2).sum()) / (A.shape[0] + A.shape[1])
errors.append(e)
print(elapsed, errors[-1])
plt.plot(errors)
plt.show()


---------------------------------------------------------------------------
ValueError                                Traceback (most recent call last)
<ipython-input-462-0ce60c95dca0> in <module>()
     26         err = a_i - U.dot(U.T)
     27         dM = - U.T.dot(err)
---> 28         M[:,m_ind:m_ind+block_size] = M[:,m_ind:m_ind+block_size] - eta * dM
     29 
     30 

ValueError: could not broadcast input array from shape (4,10) into shape (4,1)

In [ ]:


In [ ]:


In [ ]:


In [ ]:


In [431]:
U.dot(M)[3,:5]


Out[431]:
array([ 6.93616054, -3.70703094, -0.34937645, -0.53202702,  0.6691046 ])

In [457]:
a_i


Out[457]:
array([[ 0.00443189]])

In [ ]:


In [408]:
U.dot(M)[3,:5]


Out[408]:
array([ 0.15421251,  0.50652959,  0.25934971,  0.8203848 ,  0.6798117 ])

In [441]:
M


Out[441]:
array([[ 0.88991031, -1.11520329,  0.75238804,  2.13531994,  0.33131531,
         0.40815122,  0.08004733,  1.55620641,  0.0437635 ,  0.64632387,
         0.76435717, -2.19905092, -0.7758741 ,  1.16759714, -1.44449923,
         0.00943331,  1.15211436, -0.1773714 ,  0.48856296, -0.06216587],
       [ 0.0279215 ,  1.8936256 ,  0.17278202, -2.03120041, -0.83029699,
         0.24711855, -0.77822443, -0.0130944 , -0.80994036,  0.37938036,
        -1.02142506,  1.35912663, -1.00546452, -0.1222818 , -0.61429996,
         2.02352052, -2.32272918, -1.40966399,  0.6990727 , -0.24108175],
       [ 1.65758934, -1.5679221 ,  0.83967987, -0.38427077,  0.84137448,
         2.04544519,  0.02173059, -0.05216779, -0.25287257,  0.98262751,
         0.41550207, -0.17672982,  1.24510027,  0.66171195, -0.81434485,
        -0.88225348,  1.71814951, -0.37139128,  0.86622349, -0.54101977],
       [-1.92895956,  0.20296416,  0.62142512,  0.23157709,  0.43959861,
        -2.47030314,  1.03303486,  0.23278662, -0.88808655,  0.52651831,
         1.92479865, -0.64757968, -0.35915632, -0.98211636,  0.54844687,
         0.91761771,  0.58950098, -2.35283403,  0.17727504, -1.15965419]])

In [ ]: