在斯坦福大学的Coursera的Andrew Ng的机器学习介绍性讲座中，他给出了以下一行Octave解决方案的鸡尾酒会问题，因为音频源由两个空间分离的麦克风记录：

[W,s,v]=svd((repmat(sum(x.*x,1),size(x,1),1).*x)*x’);

• 传统方法

独立成分分析ICA

独立分量分析（Independent Component Analysis，ICA）是将信号之间的独立性作为分离变量判据的方法。由Comon于1994年首次提出。Comon指出ICA方法可以通过某个对比函数（Contrast Function）的目标函数达到极大值来消除观察信号中的高阶统计关联，实现盲源分离。盲源分离被描述成在不知道传输通道特性的情况下，从传感器或传感器阵列中分离或估计原始源波形的问题。然而，对于某些模型，不能保证估计或提取的信号与源信号具有完全相同的波形，因此有时要求放宽到提取的波形是源信号的失真或滤波版本。

独立成分分析（Independent Component Analysis），最早应用于盲源信号分离（Blind Source Separation，BBS）。起源于“鸡尾酒会问题”，描述如下：在嘈杂的鸡尾酒会上，许多人在同时交谈，可能还有背景音乐，但人耳却能准确而清晰的听到对方的话语。这种可以从混合声音中选择自己感兴趣的声音而忽略其他声音的现象称为“鸡尾酒会效应”。
  独立成分分析是从盲源分离技术发展而来的一种数据驱动的信号处理方法, 是基于高阶统计特性的分析方法。它利用统计原理进行计算，通过线性变换把数据或信号分离成统计独立的非高斯的信号源的线性组合。
  主成分分析（PCA）是一种数据降维的方法，它与主成分分析有一些区别.

ICA的python实现

应用Python机器学习库SKlearn中的FastICA来演示信号分离并与PCA进行对比。

#coding:utf-8
import numpy as np
import matplotlib.pyplot as plt
from scipy import signal
from sklearn.decomposition import FastICA, PCA
# 生成观测模拟数据
np.random.seed(0)
n_samples = 2000
time = np.linspace(0, 8, n_samples)
s1 = np.sin(2 * time)  # 信号源 1 : 正弦信号
s2 = np.sign(np.sin(3 * time))  # 信号源 2 : 方形信号
s3 = signal.sawtooth(2 * np.pi * time)  # 信号源 3: 锯齿波信号
S = np.c_[s1, s2, s3]
S += 0.2 * np.random.normal(size=S.shape)  # 增加噪音数据
S /= S.std(axis=0)  # 标准化

# 混合数据
A = np.array([[1, 1, 1], [0.5, 2, 1.0], [1.5, 1.0, 2.0]])  # 混合矩阵
X = np.dot(S, A.T)  # 生成观测信号源

# ICA模型
ica = FastICA(n_components=3)
S_ = ica.fit_transform(X)  # 重构信号
A_ = ica.mixing_  # 获得估计混合后的矩阵

# PCA模型
pca = PCA(n_components=3)
H = pca.fit_transform(X)  # 基于PCA的成分正交重构信号源

# 图形展示
plt.figure()
models = [X, S, S_, H]
names = ['Observations (mixed signal)',
         'True Sources',
         'ICA recovered signals',
         'PCA recovered signals']
colors = ['red', 'steelblue', 'orange']
for ii, (model, name) in enumerate(zip(models, names), 1):
    plt.subplot(4, 1, ii)
    plt.title(name)
    for sig, color in zip(model.T, colors):
        plt.plot(sig, color=color)
plt.subplots_adjust(0, 0.1, 1.2, 1.5, 0.26, 0.46)
plt.show()

from sklearn.decomposition import FastICA
ica = FastICA(n_components=None, algorithm=’parallel’, whiten=True, fun=’logcosh’, fun_args=None, max_iter=200, w_init=None)

以上是 FastICA 模型通常所含的参数及其对应的默认值。 n_components: 指定使用元素的数目。 algorithm: {‘parallel’,’deflational’},指定 FastICA 使用哪种算法。 writen: True/False，是否进行白化处理。 fun: {‘logcosh’,’exp’,’cube’,..}，选择一种近似于计算负熵的目标函数，可自己定义函数。 fun_args: 指定目标函数所要用的参数。 max_iter: 指定拟合过程中最大的迭代次数。 w_init: 指定初始的混合矩阵。

​ 1. 主成分分析假设源信号间彼此非相关，独立成分分析假设源信号间彼此独立。

​ 2. 主成分分析认为主元之间彼此正交，样本呈高斯分布；独立成分分析则不要求样本呈高斯分布。

稀疏主成分分析Spars PCA

稀疏主成分分析即是为了解决这个问题而引进的一个算法。它会把主成分系数（构成主成分时每个变量前面的系数）变的稀疏，也即是把大多数系数都变成零，通过这样一种方式，我们就可以把主成分的主要的部分凸现出来，这样主成分就会变得较为容易解释。

实现主成分分析稀疏化，最终会转化为优化问题， 也即对本来的主成分分析（PCA）中的问题增加一个惩罚函数。 这个惩罚函数包含有稀疏度的信息。当然，最终得到的问题是NP-困难问题，为了解决它，我们就需要采用一些方法逼近这个问题的解。这也是不同稀疏主成分分析算法不同的由来。

非负矩阵分解NMF

NMF的基本思想可以简单描述为：对于任意给定的一个非负矩阵V，NMF算法能够寻找到一个非负矩阵W和一个非负矩阵H，使得满足 ，从而将一个非负的矩阵分解为左右两个非负矩阵的乘积。如下图所示，其中要求分解后的矩阵H和W都必须是非负矩阵。

参考：

https://zhuanlan.zhihu.com/p/142143151

https://github.com/jake-g/dsp-fpga-labs/blob/c0c84ed08c02da11fb4161599cafe8ddcabeaca4/Labs/sound_split_demo/seperate.m

https://www.cxymm.net/article/weixin_30446557/115847250

https://leoncuhk.gitbooks.io/feature-engineering/content/feature-extracting04.html

ICA的数学推导

ICA算法的思路比较简单，但是推导过程比较复杂，本文只是梳理了推理路线。

假设我们有n个混合信号源\(X\subset{R^{n}}\)和n个独立信号\(S\subset{R^{n}}\)，且每个混合信号可以由n个独立信号的线性组合产生，即：\(X=\left[ \begin{matrix} x_1&\\ x_2&\\ …&\\ x_n& \end{matrix} \right]S=\left[ \begin{matrix} s_1&\\ s_2&\\ …&\\ s_n& \end{matrix} \right]X=AS => S=WX,W=A^{-1}\)

假设我们现在对于每个混合信号，可以取得m个样本，则有如下n*m的样本矩阵：\(D=\left[ \begin{matrix} d_{11}&d_{12}&…&d_{1m}&\\ …&\\ d_{n1}&d_{n2}&…&d_{nm}\\ \end{matrix} \right]\)

由于S中的n个独立信号是相互独立的，则它们的联合概率密度为：$$p_S(s)=\Pi_{i=1}^{n}p_{s_i}(s_i)$$

由于\(s=Wx\)，因此我们可以得出：\(p_X(x)=F^{‘}_{X}(x)=|\frac{\partial{s}}{\partial{x}}|*p_S(s(x))=|W|*\Pi_{i=1}^{n}p_{s_i}(w_ix)\)

考虑目前有m个样本，则可以得到所有样本的似然函数：\(L=\Pi_{i=1}^{m}(|W|*\Pi_{j=1}^{n}p_{s_j}(w_{j\cdot}d_{\cdot{i}}))\)

取对数之后，得到：\(lnL=\Sigma_{i=1}^{m}\Sigma_{j=1}^{n}lnp_{s_j}(w_{j\cdot}d_{\cdot{i}})+mln|W|\)

之后只要通过梯度下降法对lnL求出最大值即可，即求使得该样本出现概率最大的参数W。
此时假设我们上面的各个独立信号的概率分布函数为sigmoid函数，但是不确定这里的g函数和下面fastICA中的g函数是否有关联）：\(F_{s_i}(s_i)=\frac{1}{1+e^{-s_i}}\)

最终，我们求得：\(\frac{\partial{lnL}}{\partial{W}}=Z^TD+\frac{m}{|W|}(W^*)^T\)

其中：$$Z=g(K)=\left[ \begin{matrix} g(k_{11})&g(k_{12})&…&g(k_{1m})&\\ …&\\ g(k_{n1})&g(k_{n2})&…&g(k_{nm})\\ \end{matrix} \right]g(x)=\frac{1-e^x}{1+e^x}K=WDD=\left[ \begin{matrix} d_{11}&d_{12}&…&d_{1m}&\\ …&\\ d_{n1}&d_{n2}&…&d_{nm}\\ \end{matrix} \right]$$

由于伴随矩阵具有以下性质：$$WW^*=|W|I$$

因此我们可以求出：\(\frac{\partial{lnL}}{\partial{W}}=Z^TD+m(W^{-1})^T\)

因此可以得到梯度下降更新公式：$$W=W+\alpha(Z^TD+m(W^{-1})^T)$$

至此，ICA的基本推理就此结束。下面我们来看一下fastICA的算法过程（没有数学推理）。

fastICA的算法步骤

观测信号构成一个混合矩阵，通过数学算法进行对混合矩阵A的逆进行近似求解分为三个步骤：

• 去均值。去均值也就是中心化，实质是使信号X均值是零。
• 白化。白化就是去相关性。
• 构建正交系统。

在常用的ICA算法基础上已经有了一些改进，形成了fastICA算法。fastICA实际上是一种寻找\(w^Tz(Y=w^Tz)\)的非高斯最大的不动点迭代方案。具体步骤如下：

1. 观测数据的中心化(去均值)
2. 数据白化(去相关)，得到z
3. 选择需要顾及的独立源的个数n
4. 随机选择初始权重W（非奇异矩阵）
5. 选择非线性函数g
6. 迭代更新：
   • \(w_i \leftarrow E\{zg(w_i^Tz)\}-E\{g^{‘}(w_i^Tz)\}w\)
   • \(W \leftarrow (WW^T)^{-1/2}W\)
7. 判断收敛，是下一步，否则返回步骤6
8. 返回近似混合矩阵的逆矩阵

代码实现
基于python2.7，matplotlib，numpy实现ICA，主要参考sklean的FastICA实现。

import math
import random
import matplotlib.pyplot as plt
from numpy import *

n_components = 2

def f1(x, period = 4):
return 0.5(x-math.floor(x/period)period)

def create_data():
#data number
n = 500
#data time
T = [0.1*xi for xi in range(0, n)]
#source
S = array([[sin(xi) for xi in T], [f1(xi) for xi in T]], float32)
#mix matrix
A = array([[0.8, 0.2], [-0.3, -0.7]], float32)
return T, S, dot(A, S)

def whiten(X):
#zero mean
X_mean = X.mean(axis=-1)
X -= X_mean[:, newaxis]
#whiten
A = dot(X, X.transpose())
D , E = linalg.eig(A)
D2 = linalg.inv(array([[D[0], 0.0], [0.0, D[1]]], float32))
D2[0,0] = sqrt(D2[0,0]); D2[1,1] = sqrt(D2[1,1])
V = dot(D2, E.transpose())
return dot(V, X), V

def _logcosh(x, fun_args=None, alpha = 1):
gx = tanh(alpha * x, x); g_x = gx ** 2; g_x -= 1.; g_x *= -alpha
return gx, g_x.mean(axis=-1)

def do_decorrelation(W):
#black magic
s, u = linalg.eigh(dot(W, W.T))
return dot(dot(u * (1. / sqrt(s)), u.T), W)

def do_fastica(X):
n, m = X.shape; p = float(m); g = _logcosh
#black magic
X *= sqrt(X.shape[1])
#create w
W = ones((n,n), float32)
for i in range(n):
for j in range(i):
W[i,j] = random.random()

#compute W
maxIter = 200
for ii in range(maxIter):
    gwtx, g_wtx = g(dot(W, X))
    W1 = do_decorrelation(dot(gwtx, X.T) / p - g_wtx[:, newaxis] * W)
    lim = max( abs(abs(diag(dot(W1, W.T))) - 1) )
    W = W1
    if lim < 0.0001:
        break
return W

def show_data(T, S):
plt.plot(T, [S[0,i] for i in range(S.shape[1])], marker=”*”)
plt.plot(T, [S[1,i] for i in range(S.shape[1])], marker=”o”)
plt.show()

def main():
T, S, D = create_data()
Dwhiten, K = whiten(D)
W = do_fastica(Dwhiten)
#Sr: reconstructed source
Sr = dot(dot(W, K), D)
show_data(T, D)
show_data(T, S)
show_data(T, Sr)

if name == “main“:
main()

参考：

http://skyhigh233.com/blog/2017/04/01/ica-math/

降维——PCA

降维就是一种对高维度特征数据预处理方法。降维是将高维度的数据保留下最重要的一些特征，去除噪声和不重要的特征，从而实现提升数据处理速度的目的。在实际的生产和应用中，降维在一定的信息损失范围内，可以为我们节省大量的时间和成本。降维也成为应用非常广泛的数据预处理方法。

降维具有如下一些优点：

• 使得数据集更易使用。
• 降低算法的计算开销。
• 去除噪声。
• 使得结果容易理解。

降维的算法有很多，比如奇异值分解(SVD)、主成分分析(PCA)、因子分析(FA)、独立成分分析(ICA)。

PCA(Principal Component Analysis)，即主成分分析方法，是一种使用最广泛的数据降维算法。PCA的主要思想是将n维特征映射到k维上，这k维是全新的正交特征也被称为主成分，是在原有n维特征的基础上重新构造出来的k维特征。PCA的工作就是从原始的空间中顺序地找一组相互正交的坐标轴，新的坐标轴的选择与数据本身是密切相关的。其中，第一个新坐标轴选择是原始数据中方差最大的方向，第二个新坐标轴选取是与第一个坐标轴正交的平面中使得方差最大的，第三个轴是与第1,2个轴正交的平面中方差最大的。依次类推，可以得到n个这样的坐标轴。通过这种方式获得的新的坐标轴，我们发现，大部分方差都包含在前面k个坐标轴中，后面的坐标轴所含的方差几乎为0。于是，我们可以忽略余下的坐标轴，只保留前面k个含有绝大部分方差的坐标轴。事实上，这相当于只保留包含绝大部分方差的维度特征，而忽略包含方差几乎为0的特征维度，实现对数据特征的降维处理。

我们如何得到这些包含最大差异性的主成分方向呢？

事实上，通过计算数据矩阵的协方差矩阵，然后得到协方差矩阵的特征值特征向量，选择特征值最大(即方差最大)的k个特征所对应的特征向量组成的矩阵。这样就可以将数据矩阵转换到新的空间当中，实现数据特征的降维。

基于特征值分解协方差矩阵实现PCA算法

输入：数据集  ，需要降到k维。

1) 去平均值(即去中心化)，即每一位特征减去各自的平均值。

2) 计算协方差矩阵 ,注：这里除或不除样本数量n或n-1,其实对求出的特征向量没有影响。

3) 用特征值分解方法求协方差矩阵 的特征值与特征向量。

4) 对特征值从大到小排序，选择其中最大的k个。然后将其对应的k个特征向量分别作为行向量组成特征向量矩阵P。

5) 将数据转换到k个特征向量构建的新空间中，即Y=PX。

 PCA实例

（1）PCA的Python实现：

##Python实现PCA
import numpy as np
def pca(X,k):#k is the components you want
  #mean of each feature
  n_samples, n_features = X.shape
  mean=np.array([np.mean(X[:,i]) for i in range(n_features)])
  #normalization
  norm_X=X-mean
  #scatter matrix
  scatter_matrix=np.dot(np.transpose(norm_X),norm_X)
  #Calculate the eigenvectors and eigenvalues
  eig_val, eig_vec = np.linalg.eig(scatter_matrix)
  eig_pairs = [(np.abs(eig_val[i]), eig_vec[:,i]) for i in range(n_features)]
  # sort eig_vec based on eig_val from highest to lowest
  eig_pairs.sort(reverse=True)
  # select the top k eig_vec
  feature=np.array([ele[1] for ele in eig_pairs[:k]])
  #get new data
  data=np.dot(norm_X,np.transpose(feature))
  return data

X = np.array([[-1, 1], [-2, -1], [-3, -2], [1, 1], [2, 1], [3, 2]])

print(pca(X,1))

用sklearn的PCA与我们的PCA做个比较：

##用sklearn的PCA
from sklearn.decomposition import PCA
import numpy as np
X = np.array([[-1, 1], [-2, -1], [-3, -2], [1, 1], [2, 1], [3, 2]])
pca=PCA(n_components=1)
pca.fit(X)
print(pca.transform(X))

本以为以后再也用不到到矩阵论的知识了，没想到还是没逃过去，大四学的矩阵论现在基本上忘光了，于是又要重新拿起课本和笔记了，我感觉这次写完过段时间又会忘掉，这记性….,在学习Svd的过程中，真的是每个概念都要重新去看，比较麻烦，不说了，摸鱼去了……

奇异值分解（Singular Value Decompositionm,简称SVD）是在机器学习领域应用较为广泛的算法之一，也是学习机器学习算法绕不开的基石之一。SVD算法主要用在降维算法中的特征分解、推荐系统、自然语言处理计算机视觉等领域。奇异值分解（SVD）通俗一点讲就是将一个线性变换分解为两个线性变换，一个线性变换代表旋转，一个线性变换代表拉伸。

$$A = U\Sigma V^T \tag{1-1}$$

注：SVD是将一个矩阵分解成两个正交矩阵和一个对角矩阵，我们知道正交矩阵对应的变换是旋转变换，对角矩阵对应的变换是伸缩变换。

1、SVD算法实现

1.1 SVD原理简单回顾

有一个\(m \times n\)的实数矩阵A，我们可以将它分解成如下的形式 \( A = U\Sigma V^T \tag{1-1} \)

其中U和V均为单位正交阵，即有UU^T=I和VV^T=I，U称为左奇异矩阵，V称为右奇异矩阵，\( \Sigma\)仅在主对角线上有值，我们称它为奇异值，其它元素均为0。上面矩阵的维度分别为 $$U \in \mathbf{R}^{m\times m},\ \Sigma \in \mathbf{R}^{m\times n},\ V \in \mathbf{R}^{n\times n}$$。

正常求上面的 \( U,V,\Sigma \) 不便于求，我们可以利用如下性质

\(AA^T=U\Sigma V^TV\Sigma^TU^T=U\Sigma \Sigma^TU^T \tag{1-2}A^TA=V\Sigma^TU^TU\Sigma V^T=V\Sigma^T\Sigma V^T \tag{1-3}\)

1.2 SVD算法#

据1.1小节，对式（1-3）和式（1-4）做特征值分解，即可得到奇异值分解的结果。但是样分开求存在一定的问题，由于做特征值分解的时候，特征向量的正负号并不影响结果，比如，我们利用式（1-3）和（1-4）做特征值分解 \( AA^T\mathbf{u}_i = \sigma_i \mathbf{u}_i\quad \text{or} \quad AA^T(-\mathbf{u}_i) = \sigma_i (-\mathbf{u}_i)\\ A^TA\mathbf{v}_i = \sigma_i \mathbf{v}_i\quad \text{or} \quad A^TA(-\mathbf{v}_i) = \sigma_i (-\mathbf{v}_i)\)

如果在计算过程取，取上面的\(\mathbf{u}_i\)组成左奇异矩阵U，取 \( -\mathbf{v}_i \) 组成右奇异矩阵V，此时\(A\ne U\Sigma V^T \) 。因此求\(\mathbf{v}_i\)时，要根据\(\mathbf{u}_i\)来求，这样才能保证\(A= U\Sigma V^T\)。因此，我们可以得出如下1.1计算SVD的算法。它主要是先做特性值分解，再根据特征值分解得到的左奇异矩阵U间接地求出部分的右奇异矩阵\(V’\in \mathbf{R}^{m\times n}\)。

算法1.1：SVD输入：样本数据
输出：左奇异矩阵，奇异值矩阵，右奇异矩阵

1. 计算特征值： 特征值分解AA^T，其中\(A \in \mathbf{R}^{m\times n}\) 为原始样本数据\(AA^T=U\Sigma \Sigma^TU^T\)得到左奇异矩阵\(U \in \mathbf{R}^{m \times m}\) 和奇异值矩阵\(\Sigma’ \in \mathbf{R}^{m \times m}\)
2. 间接求部分右奇异矩阵： 求\(V’ \in \mathbf{R}^{m \times n}\)利用\(A=U\Sigma’V’\)可得\(V’ = (U\Sigma’)^{-1}A = (\Sigma’)^{-1}U^TA \tag{1-4}\)
3. 返回\(U,\ \Sigma’,\ V’\)，分别为左奇异矩阵，奇异值矩阵，右奇异矩阵。

注：这里得到的\(\Sigma’\)和V’与式（1-2）所得到的\(\Sigma,\ V\)有区别，它们的维度不一样。\(\Sigma’\)是只取了前m个奇异值形成的对角方阵，即\(\Sigma’ \in \mathbf{R}^{m \times m}；V’\)不是一个方阵，它只取了\(V \in \mathbf{R}^{m \times n}\)的前m行（假设m < n），即有\(V’ = V(:m,\cdot)\)。这样一来，我们同样有类似式（1-1）的数学关系成立，即\(A = U\Sigma’ (V’)^T\tag{1-5}\)

我们可以利用此关系重建原始数据。

2、SVD的Python实现#

以下代码的运行环境为python3.6+jupyter5.4。

2.1 SVD实现过程#

读取数据#

这里面的数据集大家随便找一个数据就好，如果有需要我的数据集，可以下在面留言。

import numpy as np
import pandas as pd
from scipy.io import loadmat

# 读取数据，使用自己数据集的路径。
train_data_mat = loadmat("../data/train_data2.mat")
train_data = train_data_mat["Data"]
print(train_data.shape)

特征值分解#

# 数据必需先转为浮点型，否则在计算的过程中会溢出，导致结果不准确
train_dataFloat = train_data / 255.0
# 计算特征值和特征向量
eval_sigma1,evec_u = np.linalg.eigh(train_dataFloat.dot(train_dataFloat.T))

计算右奇异矩阵#

#降序排列后，逆序输出
eval1_sort_idx = np.argsort(eval_sigma1)[::-1]
# 将特征值对应的特征向量也对应排好序
eval_sigma1 = np.sort(eval_sigma1)[::-1]
evec_u = evec_u[:,eval1_sort_idx]
# 计算奇异值矩阵的逆
eval_sigma1 = np.sqrt(eval_sigma1)
eval_sigma1_inv = np.linalg.inv(np.diag(eval_sigma1))
# 计算右奇异矩阵
evec_part_v = eval_sigma1_inv.dot((evec_u.T).dot(train_dataFloat))

上面的计算出的evec_u, eval_sigma1, evec_part_v分别为左奇异矩阵，所有奇异值，右奇异矩阵。

2.2 SVD降维后重建数据#

取不同个数的奇异值，重建图片，计算出均方误差，如图2-1所示。从图中可以看出，随着奇异值的增加，均方误差（MSE）在减小，且奇异值和的比率正快速上升，在100维时，奇异值占总和的53%

注： 均方误差MSE有如下计算公式\(\text{MSE} = \frac{1}{n}\left((y_1-y_1′)^2+(y_2-y_2′)^2+\cdots+(y_n-y_n’)^2\right)\).我们平时听到的\(\text{RMSE}=\sqrt{\text{MSE}}\)。

SVD与特征值分解（EVD）非常类似，应该说EVD只是SVD的一种特殊怀况。我们可以通过它们在实际的应用中返过来理解特征值/奇异值的含义：特征值/奇异值代表着数据的信息量，它的值越大，信息越多。

奇异值分解定义#

有一个\(m \times n\)的实数矩阵A，我们想要把它分解成如下的形式\(A = U\Sigma V^T \tag{2-1}\)

其中U和V均为单位正交阵，即有UU^T=I和VV^T=I，U称为左奇异矩阵，V称为右奇异矩阵，\(\Sigma\)仅在主对角线上有值，我们称它为奇异值，其它元素均为0。上面矩阵的维度分别为\(U \in R^{m\times m},\ \Sigma \in R^{m\times n},\ V \in R^{n\times n}\)。

一般地\(\Sigma\)有如下形式

\(\Sigma = \left[ \begin{matrix} \sigma_1 & 0 & 0 & 0 & 0\\ 0 & \sigma_2 & 0 & 0 & 0\\ 0 & 0 & \ddots & 0 & 0\\ 0 & 0 & 0 & \ddots & 0\\ \end{matrix} \right]_{m\times n}\)

在图像压缩中的应用

准备工具

下面的代码运行环境为python3.6+jupyter5.4

SVD（Python）

这里暂时用numpy自带的svd函数做图像压缩

①读取图片

%matplotlib inline
import matplotlib.pyplot as plt
import matplotlib.image as mpimg
import numpy as np

img_eg = mpimg.imread("../img/beauty.jpg")
print(img_eg.shape)

图片的大小是\(600\times 400 \times 3\)

②奇异值分解

img_temp = img_eg.reshape(600, 400 * 3)
U,Sigma,VT = np.linalg.svd(img_temp)

我们先将图片变成\(600\times 1200\)，再做奇异值分解。从svd函数中得到的奇异值sigma它是从大到小排列的。

③取前部分奇异值重构图片

# 取前60个奇异值
sval_nums = 60
img_restruct1 = (U[:,0:sval_nums]).dot(np.diag(Sigma[0:sval_nums])).dot(VT[0:sval_nums,:])
img_restruct1 = img_restruct1.reshape(600,400,3)

# 取前120个奇异值
sval_nums = 120
img_restruct2 = (U[:,0:sval_nums]).dot(np.diag(Sigma[0:sval_nums])).dot(VT[0:sval_nums,:])
img_restruct2 = img_restruct2.reshape(600,400,3)

将图片显示出来看一下，对比下效果

fig, ax = plt.subplots(1,3,figsize = (24,32))

ax[0].imshow(img_eg)
ax[0].set(title = "src")
ax[1].imshow(img_restruct1.astype(np.uint8))
ax[1].set(title = "nums of sigma = 60")
ax[2].imshow(img_restruct2.astype(np.uint8))
ax[2].set(title = "nums of sigma = 120")

降维——奇异值分解

奇异值分解，全称Singular Value Decomposition，简称SVD。它是一种矩阵因式分解。通过计算奇异值个数和奇异向量，生成一个可以代替原矩阵的近似矩阵，将数据集的奇异值表征按重要性排列，舍弃不重要的特征向量。可用来达到降维的目的，从而找出数据中的主成分。

奇异值分解的计算该方法类似主成分分析(PCA)，差别在于PCA利用协方差矩阵进行分解，而SVD直接在原始矩阵上进行分解。所以SVD不要求被分解的矩阵是方阵，可以操作PCA无法操作的数据集，这也是SVD有价值的特点之一。

SVD推荐系统

首先，我们知道，任何推荐系统的整体大框架都是两部分：对某个用户user而言：首先是从数百万种Item中粗略的选出千级别的Item（实现这种功能的方法叫做召回算法），然后从几百上千的Item中精细的预测算出用户对每个Item的具体评价分数，然后选出前十个评价分数最高的Item推荐给用户（实现该方法的过程称为排序算法）

基于召回的算法，我们常常使用的有两种：基于用户行为的以及基于内容的。基于用户行为的推荐算法我们也称为协同过滤算法，基于内容的推荐算法后期详细讨论。而协同过滤算法又分为两大类，一类是基于邻域的协同过滤算法，一类是基于LFM（隐因子模型），而我们今天重点讨论的SVD算法就是召回算法里的LFM算法。

我们在召回阶段的目的是由一个稀疏的评分矩阵推算出其中的空缺分数。（补充：评分矩阵是指每个用户对各个电影的评价分数，m行代表m个用户，n行代表n个电影）

我们的目的是得到用户对未评价电影的分数，而这种user-Item的评分矩阵中，用户对Item的分数是没有直接关系的，我们需要寻找一种隐空间，使得将用户和Item联系起来。比如我们无法得知user1对《利刃出鞘》这部电影的分数，但是我们可以通过用户的历史信息（即用户对别的电影的分数）得到用户对“悬疑”类电影的爱好程度，我们也可以得到某个电影的“悬疑类”的程度有多大，这样，我们就可以将这俩个关系联系到一起了，比如我们想用户A对《南方车站的故事》进行预测，我们无法直接得到他的分数，但是我们知道他对“悬疑”类的电影有多喜爱，而《南方车站的故事》的“悬疑”程度有多大，两个值相乘即可得到A对《南方车站的故事》的预测分数，这里我们只是使用到了一个维度：“悬疑”类，我们可以将矩阵A分解成m* k列的一共k个维度的矩阵，将电影分解成k * n列一共k个维度的矩阵，然后这两个矩阵相乘即可得到矩阵中空的A[i][j]的值。

说白了，我们这里使用到得SVD的思想只是说任何一个矩阵均可以分解成两个矩阵相乘的形式（SVD是三个，但也可以合并俩个得到两个矩阵），这俩个子矩阵分别代表了用户的偏好程度和电影的成分程度，将这两个矩阵相乘即可得到用户对电影的评分矩阵。

至此，我们使用SVD的思想将矩阵A分解成两个矩阵相乘的形式（B*C），但事情没有那么简单，我们的推荐系统的矩阵A是一个稀疏矩阵，基本上是不能分解的，无法直接SVD，但分解的思想我们依旧要使用，所以我们使用机器学习的方法，使用最小均方误差学习到每个矩阵的各个元素，具体来讲就是，矩阵B的i行乘矩阵C的j列得到A的ij元素，有真实值的A[i][j]元素作为标签，不断地使用RMSE减小A中所有的标签，然后得到矩阵B和C，然后使用B和C取预测A中的空白缺失值。

参考：

https://www.cnblogs.com/endlesscoding/p/10033527.html

https://www.cnblogs.com/endlesscoding/p/10058532.html

https://www.imooc.com/article/267351

我们通常使用困惑度（perplexity）来评价语言模型的好坏。困惑度是对交叉熵损失函数做指数运算后得到的值。

• 最佳情况下，模型总是把标签类别的概率预测为1，此时困惑度为1；
• 最坏情况下，模型总是把标签类别的概率预测为0，此时困惑度为正无穷；
• 基线情况下，模型总是预测所有类别的概率都相同，此时困惑度为类别个数。

困惑度（perplexity）的基本思想是：给测试集的句子赋予较高概率值的语言模型较好,当语言模型训练完之后，测试集中的句子都是正常的句子，那么训练好的模型就是在测试集上的概率越高越好

优点：

计算速度快，允许研究人员快速的淘汰不可能表现良好的模型

有助于估算 模型的不确定性和信息密度

缺点：

不适合最终评估，他只是测量模型的可信度，而不是准确性

很难在不同上下文长度、词汇大小、基于单词 与基于字符的模型等的数据集之间进行比较