embeddings = np.concatenate((vecs[my_word_indices], vecs[:10000,:]), axis=0); embeddings.shape

# (10030, 100)

from collections import defaultdict
from sklearn import decomposition

components = get_components(embeddings, categories, my_word_indices)
plotly_3d(components.T[:len(my_words),:], categories, "pca.html")

# (30, 10030)

IFrame('pca.html', width=600, height=400)

embeddings.shape

# (10030, 100)

rand_proj = embeddings @ np.random.normal(size=(embeddings.shape[1], 40)); rand_proj.shape

# (10030, 40)

# pca = decomposition.PCA(n_components=3).fit(rand_proj.T)
# components = pca.components_
components = get_components(rand_proj, categories, my_word_indices)
plotly_3d(components.T[:len(my_words),:], categories, "pca-rand-proj.html")

# (30, 10030)

IFrame('pca-rand-proj.html', width=600, height=400)

import imageio
imageio.plugins.ffmpeg.download()

import moviepy.editor as mpe
import numpy as np
import scipy

%matplotlib inline
import matplotlib.pyplot as plt

scale = 0.50   # adjust the ratio to change the resolution of the image
dims = (int(240 * scale), int(320 * scale))
fps = 60      # frames per second

M = np.load("movie/med_res_surveillance_matrix_60fps.npy")

print(dims, M.shape)

# (120, 160) (19200, 6000)

plt.imshow(np.reshape(M[:,140], dims), cmap='gray');

plt.figure(figsize=(12, 12))
plt.imshow(M, cmap='gray')

# <matplotlib.image.AxesImage at 0x7f601f315fd0>

U = np.load("U.npy")
s = np.load("s.npy")
V = np.load("V.npy")

U.shape, s.shape, V.shape

# ((19200, 6000), (6000,), (6000, 6000))

reconstructed_matrix = U @ np.diag(s) @ V

np.allclose(M, reconstructed_matrix)

# True

low_rank = np.expand_dims(U[:,0], 1) * s[0] * np.expand_dims(V[0,:], 0)

plt.figure(figsize=(12, 12))
plt.imshow(low_rank, cmap='gray')

# <matplotlib.image.AxesImage at 0x7f1cc3e2c9e8>

plt.imshow(np.reshape(low_rank[:,0], dims), cmap='gray');

plt.imshow(np.reshape(M[:,0] - low_rank[:,0], dims), cmap='gray');

np.set_printoptions(suppress=True, precision=4)

import timeit
import pandas as pd

m_array = np.array([100, int(1e3), int(1e4)])
n_array = np.array([100, int(1e3), int(1e4)])

index = pd.MultiIndex.from_product([m_array, n_array], names=['# rows', '# cols'])

pd.options.display.float_format = '{:,.3f}'.format
df = pd.DataFrame(index=m_array, columns=n_array)

# %%prun
for m in m_array:
    for n in n_array:      
        A = np.random.uniform(-40,40,[m,n])  
        t = timeit.timeit('np.linalg.svd(A, full_matrices=False)', number=3, globals=globals())
        df.set_value(m, n, t)

df/3

%time u, s, v = np.linalg.svd(M, full_matrices=False)

'''
CPU times: user 5min 38s, sys: 1.53 s, total: 5min 40s
Wall time: 57.1 s
'''

M.shape

# (19200, 6000)

def simple_randomized_svd(M, k=10):
    m, n = M.shape
    transpose = False
    if m < n:
        transpose = True
        M = M.T
        
    rand_matrix = np.random.normal(size=(M.shape[1], k))  # short side by k
    Q, _ = np.linalg.qr(M @ rand_matrix, mode='reduced')  # long side by k
    smaller_matrix = Q.T @ M                              # k by short side
    U_hat, s, V = np.linalg.svd(smaller_matrix, full_matrices=False)
U = Q @ U_hat

if transpose:
    return V.T, s.T, U.T
else:
    return U, s, V

%time u, s, v = simple_randomized_svd(M, 10)

**CPU times: user 3.06 s, sys: 268 ms, total: 3.33 s
Wall time: 789 ms**

U_rand, s_rand, V_rand = simple_randomized_svd(M, 10)

low_rank = np.expand_dims(U_rand[:,0], 1) * s_rand[0] * np.expand_dims(V_rand[0,:], 0)

plt.imshow(np.reshape(low_rank[:,0], dims), cmap='gray');

Мы как получаем людей?

Как этот метод работает?
```py
rand_matrix = np.random.normal(size=(M.shape[1], 10))

rand_matrix.shape
# (6000, 10)

plt.imshow(np.reshape(rand_matrix[:4900,0], (70,70)), cmap='gray');

temp = M @ rand_matrix; temp.shape
# (19200, 10)

plt.imshow(np.reshape(temp[:,0], dims), cmap='gray');

plt.imshow(np.reshape(temp[:,1], dims), cmap='gray');

Q, _ = np.linalg.qr(M @ rand_matrix, mode='reduced'); Q.shape
# (19200, 10)

np.dot(Q[:,0], Q[:,1])
# -3.8163916471489756e-17

plt.imshow(np.reshape(Q[:,0], dims), cmap='gray');

plt.imshow(np.reshape(Q[:,1], dims), cmap='gray');

smaller_matrix = Q.T @ M; smaller_matrix.shape
# (10, 6000)

U_hat, s, V = np.linalg.svd(smaller_matrix, full_matrices=False)

U = Q @ U_hat

plt.imshow(np.reshape(U[:,0], dims), cmap='gray');

reconstructed_small_M = U @ np.diag(s) @ V

plt.imshow(np.reshape(M[:,0] - reconstructed_small_M[:,0], dims), cmap='gray');

from sklearn import decomposition
import fbpca

%time u, s, v = np.linalg.svd(M, full_matrices=False)
'''
CPU times: user 5min 38s, sys: 1.53 s, total: 5min 40s
Wall time: 57.1 s
'''

%time u, s, v = simple_randomized_svd(M, 10)
'''
CPU times: user 2.37 s, sys: 160 ms, total: 2.53 s
Wall time: 641 ms
'''

%time u, s, v = decomposition.randomized_svd(M, 10)
'''
CPU times: user 19.2 s, sys: 1.44 s, total: 20.7 s
Wall time: 3.67 s
'''

import timeit
import pandas as pd

U_rand, s_rand, V_rand = fbpca.pca(M, 700, raw=True)
reconstructed = U_rand @ np.diag(s_rand) @ V_rand

np.linalg.norm(M - reconstructed)
# 1.1065914828881536e-07

plt.imshow(np.reshape(reconstructed[:,140], dims), cmap='gray');

pd.options.display.float_format = '{:,.2f}'.format
k_values = np.arange(100,1000,100)
df_rand = pd.DataFrame(index=["time", "error"], columns=k_values)

# df_rand = pd.read_pickle("svd_df")

for k in k_values:
    U_rand, s_rand, V_rand = fbpca.pca(M, k, raw=True)
    reconstructed = U_rand @ np.diag(s_rand) @ V_rand
    df_rand.set_value("error", k, np.linalg.norm(M - reconstructed))
    t = timeit.timeit('fbpca.pca(M, k)', number=3, globals=globals())
    df_rand.set_value("time", k, t/3)

df_rand.to_pickle("df_rand")

df_rand

df = pd.DataFrame(index=["error"], columns=k_values)

for k in k_values:
    reconstructed = U[:,:k] @ np.diag(s[:k]) @ V[:k,:]
    df.set_value("error", k, np.linalg.norm(M - reconstructed))

df.to_pickle("df")

fig, ax1 = plt.subplots()
ax1.plot(df.columns, df_rand.loc["time"].values, 'b-', label="randomized SVD time")
ax1.plot(df.columns, np.tile([57], 9), 'g-', label="SVD time")
ax1.set_xlabel('k: # of singular values')
# Make the y-axis label, ticks

В этом тексте описывается процесс работы с матрицами и их разложениями с использованием различных методов и библиотек. Текст содержит фрагменты кода на языке Python, а также результаты вычислений и сравнения времени выполнения различных алгоритмов. **Метод randomized_range_finder находит ортогональную матрицу, диапазон которой приблизительно равен диапазону матрицы A (шаг 1 нашего алгоритма).**

Для этого мы используем LU и QR разложения, которые рассмотрим позже.

Я использовал исходный код sklearn.extmath.randomized_svd в качестве руководства.
```py
# Вычисление ортогональной матрицы, диапазон которой примерно равен диапазону A
# power_iteration_normalizer может быть safe_sparse_dot (быстрый, но нестабильный), LU (средний) или QR (медленный, но наиболее точный)
def randomized_range_finder(A, size, n_iter=5):
    Q = np.random.normal(size=(A.shape[1], size))

    for i in range(n_iter):
        Q, _ = linalg.lu(A @ Q, permute_l=True)
        Q, _ = linalg.lu(A.T @ Q, permute_l=True)

    Q, _ = linalg.qr(A @ Q, mode='economic')
    return Q

def randomized_svd(M, n_components, n_oversamples=10, n_iter=4):

    n_random = n_components + n_oversamples

    Q = randomized_range_finder(M, n_random, n_iter)
    print(Q.shape)
    # проецируем M на (k + p)-мерное пространство с использованием базисных векторов
    B = Q.T @ M
    print(B.shape)
    # вычисляем SVD на тонкой матрице: ширина (k + p)
    Uhat, s, V = linalg.svd(B, full_matrices=False)
    del B
    U = Q @ Uhat
    print(U.shape)

    return U[:, :n_components], s[:n_components], V[:n_components, :]

u, s, v = randomized_svd(vectors, 5)

vectors.shape

# (2034, 26576)

Q = np.random.normal(size=(vectors.shape[1], 10)); Q.shape

# (26576, 10)

Q2, _ = linalg.qr(vectors @ Q, mode='economic'); Q2.shape

# (2034, 10)

Q2.shape

# (2034, 10)

%time u, s, v = randomized_svd(vectors, 5)

'''
CPU times: user 136 ms, sys: 0 ns, total: 136 ms
Wall time: 137 ms
'''

u.shape, s.shape, v.shape

# ((2034, 5), (5,), (5, 26576))

show_topics(v)

'''
['jpeg image edu file graphics images gif data',
 'edu graphics data space pub mail 128 3d',
 'space jesus launch god people satellite matthew atheists',
 'space launch satellite commercial nasa satellites market year',
 'image data processing analysis software available tools display']
'''

plt.plot(range(0,n*step,step), error)

# [<matplotlib.lines.Line2D at 0x7fe3f8a1b438>]

%time u, s, v = decomposition.randomized_svd(vectors, 5)

'''
CPU times: user 144 ms, sys: 8 ms, total: 152 ms
Wall time: 154 ms
'''

%time u, s, v = decomposition.randomized_svd(vectors.todense(), 5)

'''
CPU times: user 2.38 s, sys: 592 ms, total: 2.97 s
Wall time: 2.96 s
'''