<SYSTEM_TASK:>
Given the following text description, write Python code to implement the functionality described below step by step
<END_TASK>
<USER_TASK:>
Description:
 Matrix Factorization with SimEc
In this notebook are some examples to show that SimEc can be used to compute the SVD or eigendecomposition of a matrix.
Performing an SVD and eigendecomposition of a matrix with neural networks was first described by A. Cichocki (et al.) in 1992 in this paper (SVD) and this paper (eigendecomposition).
Step1:  SVD of a random matrix
Step2:  Dealing with missing values
A regular SVD can not be computed for a matrix with missing values, therefore these need to be filled by some values, e.g. the average of the matrix. SimEc on the other hand can compute the backpropagation error considering only the available entries of the matrix, which means it can easily handle about 50% missing values without seeing major performance decreases and consistently outperforms the regular SVD (results would probably be even better with more careful hyperparameter tuning, e.g. using some regularization).
Step3:  Eigendecomposition of a square symmetric matrix
<ASSISTANT_TASK:>
Python Code:
from __future__ import unicode_literals, division, print_function, absolute_import
import numpy as np
np.random.seed(28)
import matplotlib.pyplot as plt
import tensorflow as tf
tf.set_random_seed(28)
import keras
from simec import SimilarityEncoder
%matplotlib inline
%load_ext autoreload
%autoreload 2
def msqe(A, B):
    # compute the mean squared error between two matrices A and B
    # obviously, A and B should be the same shape...
    return np.mean((A - B) ** 2)
# generate a random matrix
n_input = 500
n_output = 700
A = np.random.rand(n_input, n_output)
# compute its SVD
U, s, Vh = np.linalg.svd(A, full_matrices=True)
U.shape, s.shape, Vh.shape
# make the eigenvalues of A a bit more extreme
S = np.zeros((n_input, n_output))
s[0] = s[1]+1
s[:10] *= 50.
s[10:20] *= 20.
s[20:100] *= 10.
S[:n_input, :n_input] = np.diag(s)
# recompute A and scale it to be in a somewhat reasonable range
A = np.dot(U, np.dot(S, Vh))
A = A/np.max(np.abs(A))
# recompute SVD again
U, s, Vh = np.linalg.svd(A, full_matrices=True)
S = np.zeros((n_input, n_output))
S[:n_input, :n_input] = np.diag(s)
# check that eigenvectors are orthogonal
np.dot(Vh[:100,:], Vh[:100,:].T)
# inspect eigenvalues
plt.plot(s)
Explanation: Matrix Factorization with SimEc
In this notebook are some examples to show that SimEc can be used to compute the SVD or eigendecomposition of a matrix.
Performing an SVD and eigendecomposition of a matrix with neural networks was first described by A. Cichocki (et al.) in 1992 in this paper (SVD) and this paper (eigendecomposition).
End of explanation
# mean squared error of approximation decreases with more embedding dim
for e_dim in [2, 10, 25, 50, 75, 100, 250, 400, 500]:
    print("mse with %3i e_dim: %.8f" % (e_dim, msqe(A, np.dot(U[:,:e_dim], np.dot(S[:e_dim,:e_dim], Vh[:e_dim,:])))))
# factorize the matrix with a simec
X = np.eye(n_input)
mses = []
e_dims =  [2, 10, 25, 50, 75, 100, 250, 400, 500, 750, 1000]
l_rates = {400: 0.004, 500: 0.0034, 750: 0.0032, 1000: 0.003}
for e_dim in e_dims:
    model = SimilarityEncoder(n_input, e_dim, n_output, orth_reg=0.001 if e_dim > 500 else 0., opt=keras.optimizers.Adamax(lr=0.005 if e_dim <= 250 else l_rates[e_dim]))
    model.fit(X, A, epochs=50)
    mse = msqe(A, model.predict(X))
    mses.append(mse)
    print("mse with %4i e_dim: %.8f" % (e_dim, mse))
for i, e_dim in enumerate(e_dims):
    print("mse with %4i e_dim: %.8f" % (e_dim, mses[i]))
# factorize the matrix with a simec - transpose works just as well
# this time we use the embedding we learned before as the weights of the last layer
# to get the mapping function for the other side of the equation (which is of course kind of useless
# here since we don't actually map feature vectors but only identiy vectors)
X1 = np.eye(n_input)
X2 = np.eye(n_output)
mses1 = []
mses2 = []
mses3 = []
e_dims =  [2, 10, 25, 50, 75, 100, 250, 400, 500, 750, 1000]
for e_dim in e_dims:
    model = SimilarityEncoder(n_input, e_dim, n_output, opt=keras.optimizers.Adamax(lr=0.005 if e_dim <= 250 else l_rates[e_dim]))
    model.fit(X1, A, epochs=50)
    mse = msqe(A, model.predict(X1))
    mses1.append(mse)
    print("mse with %4i e_dim: %.8f" % (e_dim, mse))
    Y1 = model.transform(X1)
    model = SimilarityEncoder(n_output, e_dim, n_input, W_ll=Y1.T, wll_frozen=True, opt=keras.optimizers.Adamax(lr=0.005 if e_dim <= 250 else l_rates[e_dim]))
    model.fit(X2, A.T, epochs=50)
    mse = msqe(A.T, model.predict(X2))
    mses2.append(mse)
    print("mse with %4i e_dim: %.8f" % (e_dim, mse))
    Y2 = model.transform(X2)
    # the dot product of both embeddings should also approximate A
    mse = msqe(A, np.dot(Y1, Y2.T))
    mses3.append(mse)
for i, e_dim in enumerate(e_dims):
    print("mse with %4i e_dim: %.8f / %.8f / %.8f" % (e_dim, mses1[i], mses2[i], mses3[i]))
Explanation: SVD of a random matrix
End of explanation
X = np.eye(n_input)
missing_targets = [0., 0.1, 0.2, 0.3, 0.4, 0.5, 0.6, 0.7, 0.8, 0.9]
e_dims =  [2, 10, 50, 100, 250, 500]
for e_dim in e_dims:
    np.random.seed(15)
    mses = []
    mses_svd = []
    mse_svd = msqe(A, np.dot(U[:,:e_dim], np.dot(S[:e_dim,:e_dim], Vh[:e_dim,:])))
    for m in missing_targets:
        print(m)
        A_noisy = A.copy()
        A_noisy[np.random.rand(*A_noisy.shape)<=m] = -100
        model = SimilarityEncoder(n_input, e_dim, n_output, mask_value=-100, l2_reg_emb=0.00001, 
                                  l2_reg_out=0. if m < 0.7 else 0.00001, opt=keras.optimizers.Adamax(lr=0.025 if e_dim < 50 else 0.01))
        model.fit(X, A_noisy, epochs=60)
        mse = msqe(A, model.predict(X))
        mses.append(mse)
        A_noisy[A_noisy == -100] = np.mean(A)
        U_n, s_n, Vh_n = np.linalg.svd(A_noisy, full_matrices=True)
        S_n = np.zeros((n_input, n_output))
        S_n[:n_input, :n_input] = np.diag(s_n)
        mses_svd.append(msqe(A, np.dot(U_n[:,:e_dim], np.dot(S_n[:e_dim,:e_dim], Vh_n[:e_dim,:]))))
    print(mses)
    plt.figure();
    plt.plot([0, missing_targets[-1]], [mse_svd, mse_svd], '--', linewidth=0.5, label='SVD noise free');
    plt.plot(missing_targets, mses_svd, '-o', markersize=3, label='SVD');
    plt.plot(missing_targets, mses, '-o', markersize=3, label='SimEc SVD');
    plt.legend(loc=0);
    plt.title('Matrix factorization of A (%i embedding dim)' % e_dim);
    plt.xticks(missing_targets, missing_targets);
    plt.xlabel('Fraction of Missing Entries');
    plt.ylabel('Mean Squared Error');
Explanation: Dealing with missing values
A regular SVD can not be computed for a matrix with missing values, therefore these need to be filled by some values, e.g. the average of the matrix. SimEc on the other hand can compute the backpropagation error considering only the available entries of the matrix, which means it can easily handle about 50% missing values without seeing major performance decreases and consistently outperforms the regular SVD (results would probably be even better with more careful hyperparameter tuning, e.g. using some regularization).
End of explanation
# get 2 square symmetric matrices as AA^T and A^TA
S1 = np.dot(A, A.T)
S2 = np.dot(A.T, A)
# make sure their range of values is still sort of reasonable (for >100 we should probably rescale)
np.max(S1), np.min(S1), np.max(S2), np.min(S2)
# the corresponding SVD eigenvalues and -vectors should work here as well
for e_dim in [2, 10, 25, 50, 75, 100, 250, 400, 500]:
    print("mse with %3i e_dim: %11.8f" % (e_dim, msqe(S1, np.dot(U[:,:e_dim], np.dot(S[:e_dim,:e_dim]**2, U.T[:e_dim,:])))))
S_out = np.zeros((n_output, n_output))
S_out[:n_input, :n_input] = S[:n_input, :n_input]
for e_dim in [2, 10, 25, 50, 75, 100, 250, 500, 700]:
    print("mse with %3i e_dim: %11.8f" % (e_dim, msqe(S2, np.dot(Vh.T[:,:e_dim], np.dot(S_out[:e_dim,:e_dim]**2, Vh[:e_dim,:])))))
# factorize the similarity matrix S1 with a simec
X = np.eye(n_input)
mses1 = []
mses2 = []
e_dims =  [2,   10,  25,  50,  75,  100,  250,  400,  500,  750,  1000]
l_rates = [0.1, 0.1, 0.1, 0.1, 0.1, 0.05, 0.05, 0.01, 0.01, 0.01, 0.01]
for i, e_dim in enumerate(e_dims):
    model = SimilarityEncoder(n_input, e_dim, S1.shape[1], s_ll_reg=1., S_ll=S1, opt=keras.optimizers.Adamax(lr=l_rates[i]))
    model.fit(X, S1, epochs=100)
    mse = msqe(S1, model.predict(X))
    mses1.append(mse)
    print("mse with %4i e_dim: %11.8f" % (e_dim, mse))
    # scalar product of emedding should also approximate S1
    Y = model.transform(X)
    mse = msqe(S1, np.dot(Y, Y.T))
    mses2.append(mse)
    print("mse with %4i e_dim: %11.8f" % (e_dim, mse))
for i, e_dim in enumerate(e_dims):
    print("mse with %4i e_dim: %11.8f (%11.8f)" % (e_dim, mses1[i], mses2[i]))
# factorize the similarity matrix S2 with a simec
X = np.eye(n_output)
mses1 = []
mses2 = []
e_dims =  [2,   10,  25,  50,  75,  100,  250,  400,  500,  750,  1000]
l_rates = [0.1, 0.1, 0.1, 0.1, 0.1, 0.05, 0.05, 0.01, 0.01, 0.01, 0.01]
for i, e_dim in enumerate(e_dims):
    model = SimilarityEncoder(n_output, e_dim, S2.shape[1], s_ll_reg=1., S_ll=S2, opt=keras.optimizers.Adamax(lr=l_rates[i]))
    model.fit(X, S2, epochs=60)
    mse = msqe(S2, model.predict(X))
    mses1.append(mse)
    print("mse with %4i e_dim: %11.8f" % (e_dim, mse))
    # scalar product of emedding should also approximate S1
    Y = model.transform(X)
    mse = msqe(S2, np.dot(Y, Y.T))
    mses2.append(mse)
    print("mse with %4i e_dim: %11.8f" % (e_dim, mse))
for i, e_dim in enumerate(e_dims):
    print("mse with %4i e_dim: %11.8f (%11.8f)" % (e_dim, mses1[i], mses2[i]))
Explanation: Eigendecomposition of a square symmetric matrix
End of explanation
<END_TASK>