一 SVD的數學模型
設有樣本資料 X = ( x i j ) n ∗ p ; i = 1... p , j = 1... n , ∀ x i j ∈ R X=(x_{ij})_{n*p}; i=1...p,j=1...n ,\forall{x_{ij}}\in \mathbb R X=(xij)n∗p;i=1...p,j=1...n,∀xij∈R ; n n n為樣本數, p p p是樣本的觀察名額數目。将樣本資料中的一行也就是一個觀察樣本看成是P維空間的一個點, n n n個樣本也就是P維空間中的 n n n點.
矩陣 X X X有如下分解: X = U Σ V t X=U\Sigma V^t X=UΣVt,其中 U ∈ R n ∗ n U \in R^{n*n} U∈Rn∗n, V ∈ R p ∗ p V\in R^{p*p} V∈Rp∗p,$\Sigma \in R^{n*p},v_{ij}=0 , \forall{ i\neq j } , 并 且 ,并且 ,并且S, D$均為正交陣。
1.1 SVD與PCA的關系
事實上,隻要找到 H = X t X H=X^{t}X H=XtX的特征向量矩陣V,令 Y = X ∗ V Y=X*V Y=X∗V,則有 Y t Y = Λ Y^tY=\Lambda YtY=Λ, Λ \Lambda Λ的非對角線元素全為零,即Y個各列之間正交,隻需要将Y的每列進行機關化就可以得到機關正交陣的p個列,補齊其餘n-p列得到機關正交陣 U U U,即 Y = U Σ Y=U\Sigma Y=UΣ,是以有 X = U Σ V t X=U\Sigma V^t X=UΣVt成立。
首先,令 H = X t ∗ X H=X^t*X H=Xt∗X,根據SVD分解,有 H = U Σ t V t V Σ U t = U Σ t Σ U t = U t Λ U H=U\Sigma^t V^tV\Sigma U^t=U\Sigma^t \Sigma U^t=U^t\Lambda U H=UΣtVtVΣUt=UΣtΣUt=UtΛU。可見,U是H的特征向量矩陣,特征值為 λ i = σ i i 2 \lambda _{i}=\sigma_{ii}^2 λi=σii2。
其次,令 H ∗ = X ∗ X t H^*=X*X^t H∗=X∗Xt,根據SVD分解,有 H = V Σ Σ t V t = V t Λ ∗ V H=V\Sigma \Sigma^t V^t=V^t\Lambda^* V H=VΣΣtVt=VtΛ∗V。可見,V是 H ∗ H^* H∗的特征向量矩陣,特征值根為 λ i ∗ = σ i i 2 \lambda _{i}^*=\sigma_{ii}^2 λi∗=σii2。
由此知道任意的實數矩陣 X X X, 對稱陣 X ∗ X t , X t ∗ X X*X^t,X^t*X X∗Xt,Xt∗X有相同的非零特征值,他們的特征向量矩陣,對應于SVD分解的U和V。
1.2 最優化方法看SVD
用 u ⃗ , v ⃗ , λ \vec u ,\vec v,\lambda u
,v
,λ逼近 X X X,建立如下最優化方法:
m i n ( Q ( X , u ⃗ , v ⃗ , λ ) ) = m i n ( ∑ i = 1 n ∑ j = 1 p ( x i j − λ u i v j ) 2 ) s t . { ∑ i = 1 n u i 2 = 1 ∑ j = 1 p v i 2 = 1 min( Q(X,\vec u ,\vec v,\lambda))=min( \sum_{i=1}^n \sum_{j=1}^p(x_{ij}-\lambda u_i v_j)^2) \\ st. \begin{cases} \sum_{i=1}^n u_i^2 &=1\\ \\ \sum_{j=1}^p v_i^2 & =1 \end{cases} min(Q(X,u
,v
,λ))=min(i=1∑nj=1∑p(xij−λuivj)2)st.⎩⎪⎨⎪⎧∑i=1nui2∑j=1pvi2=1=1
運用拉格朗日方法,可以得到下面的方程組:
{ u ⃗ t X v ⃗ = λ X v ⃗ = λ u ⃗ ( 1 ) X t u ⃗ = λ v ⃗ \begin{cases} \vec u^t X \vec v & =\lambda\\ X \vec v & =\lambda \vec u & & & & (1) \\ X ^t\vec u &=\lambda \vec v \end{cases} ⎩⎪⎨⎪⎧u
tXv
Xv
Xtu
=λ=λu
=λv
(1)
由(1)可以得到,$XX ^t\vec u =\lambda X\vec v =\lambda^2 \vec u KaTeX parse error: Expected 'EOF', got '&' at position 2: ,&̲nbsp; 容易看到 \vec u 是 是 是XX ^t$的特征向量。
假設有兩組解滿足(1)式,分别為 ( u ⃗ 1 , v ⃗ 1 , λ 1 ) , ( u ⃗ 2 , v ⃗ 2 , λ 2 ) (\vec u _1,\vec v_1,\lambda_1),(\vec u _2,\vec v_2,\lambda_2) (u
1,v
1,λ1),(u
2,v
2,λ2),且有 λ 1 2 ≠ λ 2 2 \lambda_1^2 \neq \lambda_2^2 λ12̸=λ22則有:
u ⃗ 1 t X X t u ⃗ 2 = λ 1 2 u ⃗ 1 t u ⃗ 2 = λ 2 2 u ⃗ 1 t u ⃗ 2 \vec u _1^t XX ^t \vec u _2=\lambda_1^2 \vec u _1^t\vec u _2= \lambda_2^2 \vec u _1^t\vec u _2 u
1tXXtu
2=λ12u
1tu
2=λ22u
1tu
2,則必有 u ⃗ 1 t u ⃗ 2 = 0 \vec u _1^t\vec u _2=0 u
1tu
2=0,同理還有 v ⃗ 1 t v ⃗ 2 = 0 \vec v_1^t\vec v_2=0 v
1tv
2=0
滿足(1)式的所有解,構成了X的SVD分解。
1.3 SVD分解算法
事實上(1)給出了一個疊代算法求解SVD;基本步驟如下:
-
初始化 v ⃗ = ( 1 , 1 , … , 1 ) t \vec v=(1,1,…,1)^t v
=(1,1,…,1)t, λ 1 = 0 \lambda_1=0 λ1=0,指定一個比較小的數值 e , e > 0 e,e>0 e,e>0
-
計算 X v ⃗ = u ⃗ ′ X \vec v =\vec u' Xv
=u
′, u ⃗ = u ⃗ ′ / ∣ u ⃗ ′ ∣ \vec u ={\vec u'}/{|\vec u'|} u
=u
′/∣u
′∣
-
計算 v ⃗ ′ = X t u ⃗ \vec v'=X^t\vec u v
′=Xtu
, λ 2 = ∣ v ⃗ ′ ∣ \lambda_2={|\vec v'|} λ2=∣v
′∣, v ⃗ = v ⃗ ′ / λ 2 \vec v ={\vec v'}/\lambda_2 v
=v
′/λ2
-
計算 判斷 ∣ λ 2 − λ 1 ∣ < e |\lambda_2-\lambda_1|<e ∣λ2−λ1∣<e,為真,傳回 ( λ 2 , u ⃗ , v ⃗ ) (\lambda_2,\vec u,\vec v) (λ2,u
,v
),否則 λ 1 ← λ 2 \lambda_1 \leftarrow \lambda_2 λ1←λ2跳轉到第2步。
上面給出的分解SVD算法,隻給出了一組解,為了求得剩餘的解,隻需要将X做如下的變換: X : = X − λ ∗ u ⃗ ∗ v ⃗ t X:=X-\lambda * \vec u *\vec v^t X:=X−λ∗u
∗v
t然後按照上面的算法再計算即可,
#####下面驗證第二組解滿足(1)式:
事實上,假設兩輪計算得到了兩組解: ( λ 1 , u ⃗ 1 , v ⃗ 1 ) , ( λ 2 , u ⃗ 2 , v ⃗ 2 ) (\lambda_1,\vec u_1,\vec v_1),(\lambda_2,\vec u_2,\vec v_2) (λ1,u
1,v
1),(λ2,u
2,v
2),先驗證若: λ 1 2 ≠ λ 2 2 \lambda_1^2 \neq \lambda_2^2 λ12̸=λ22有 u ⃗ 1 , u ⃗ 2 \vec u_1,\vec u_2 u
1,u
2正交。
( X t − λ 1 ∗ v ⃗ 1 ∗ u ⃗ 1 t ) u ⃗ 1 = X t u ⃗ 1 − λ 1 ∗ v ⃗ 1 = 0 (X^t-\lambda_1 * \vec v_1*\vec u_1^t)\vec u_1=X^t\vec u_1-\lambda_1 * \vec v_1=0 (Xt−λ1∗v
1∗u
1t)u
1=Xtu
1−λ1∗v
1=0
是以有:
0 = u ⃗ 2 t ( X − λ 1 ∗ u ⃗ 1 ∗ v ⃗ 1 t ) ( X t − λ 1 ∗ v ⃗ 1 ∗ u ⃗ 1 t ) u ⃗ 1 = λ 2 2 u ⃗ 2 t u ⃗ 1 0=\vec u_2^t(X-\lambda_1 * \vec u_1 *\vec v_1^t)(X^t-\lambda_1 * \vec v_1*\vec u_1^t)\vec u_1= \lambda_2^2\vec u_2^t\vec u_1 0=u
2t(X−λ1∗u
1∗v
1t)(Xt−λ1∗v
1∗u
1t)u
1=λ22u
2tu
1
可見: u ⃗ 2 t u ⃗ 1 = 0 \vec u_2^t\vec u_1=0 u
2tu
1=0,同樣,還有: v ⃗ 2 t v ⃗ 1 = 0 \vec v_2^t\vec v_1=0 v
2tv
1=0,是以有: λ 2 = u ⃗ 2 t ( X − λ 1 ∗ u ⃗ 1 ∗ v ⃗ 1 t ) v ⃗ 2 = u ⃗ 2 t X v ⃗ 2 − λ 1 u ⃗ 2 t ∗ u ⃗ 1 ∗ v ⃗ 1 t v ⃗ 2 = u ⃗ 2 t X v ⃗ 2 λ 2 v ⃗ 2 = ( X − λ 1 ∗ u ⃗ 1 ∗ v ⃗ 1 t ) t u ⃗ 2 = X t u ⃗ 2 λ 2 u ⃗ 2 = ( X − λ 1 ∗ u ⃗ 1 ∗ v ⃗ 1 t ) v ⃗ 2 = X v ⃗ 2 \lambda_2 =\vec u_2^t(X-\lambda_1 * \vec u_1 *\vec v_1^t)\vec v_2=\vec u_2^tX\vec v_2-\lambda_1\vec u_2^t * \vec u_1 *\vec v_1^t\vec v_2=\vec u_2^tX\vec v_2 \\ \lambda_2 \vec v_2=(X-\lambda_1 * \vec u_1 *\vec v_1^t)^t\vec u_2=X^t\vec u_2\\ \lambda_2 \vec u_2=(X-\lambda_1 * \vec u_1 *\vec v_1^t)\vec v_2=X\vec v_2 λ2=u
2t(X−λ1∗u
1∗v
1t)v
2=u
2tXv
2−λ1u
2t∗u
1∗v
1tv
2=u
2tXv
2λ2v
2=(X−λ1∗u
1∗v
1t)tu
2=Xtu
2λ2u
2=(X−λ1∗u
1∗v
1t)v
2=Xv
2
可見: ( λ 1 , u ⃗ 1 , v ⃗ 1 ) , ( λ 2 , u ⃗ 2 , v ⃗ 2 ) (\lambda_1,\vec u_1,\vec v_1),(\lambda_2,\vec u_2,\vec v_2) (λ1,u
1,v
1),(λ2,u
2,v
2)都是X的SVD分解的特征。
事實上,可以通過計算 H = X t ∗ X H=X^t*X H=Xt∗X,然後用上面的辦法求解H的特征分解,算法可以進一步優化。但是當H或者X資料量特别大的時候,需要進行并行加速。上述算法可以進行并行運算,隻需要對X進行分塊即可。
二 、SVD的應用
SVD分解的應用十分廣泛,在資料壓縮,推薦,經濟結構分析等等方面,并且在應用方面與因子分析、對應分析方法極為相似。
2.1 資料壓縮
在大資料時代,存在大型數值矩陣,比如客戶與購買的商品關系矩陣,網名點選網站資料等等,樣本資料量為 n ∗ p n*p n∗p,而樣本矩陣的秩 R ( X ) = r R(X)=r R(X)=r,則SVD分解之後的去掉零後隻剩 ( n + p + 1 ) ∗ r (n+p+1)*r (n+p+1)∗r個資料需要存儲。這是還是無損壓縮。
假設 n = 100 ∗ 1 0 5 , p = 1 0 5 , r = 1 0 3 n=100*10^5,p=10^5,r=10^3 n=100∗105,p=105,r=103,則存儲原始資料量為 1 0 12 10^{12} 1012,而通過SVD分解後,存儲資料量為 ( 100 ∗ 1 0 5 + 1 0 5 + 1 ) ∗ 1 0 3 ≈ 1 0 10 (100*10^5+10^5+1)*10^3 \approx10^{10} (100∗105+105+1)∗103≈1010,所需存儲空間縮小兩個數量級。
除此之外,參照PCA(主成分分析),可将特征值從大到小排列,取方差和占比前80%的特征分量,還可以進一步壓縮資料,降低噪音。
2.2 缺失值估計,商品推薦
在很多時候,X是稀疏矩陣,比如若X的每行代表某人購買不同商品數量或者金額的樣本資料。将樣本中缺失的資料視為潛在的消費需求,那麼估計樣本中缺失資料的值,就可以用來做商品推薦。
1.2節中讨論了二乘法逼近樣本的最優化方法,并且從中給出了SVD分解的算法,現在稍微修改一下該方法,在此之前定義一個函數:
s ( x ) = { 0 x = 0 1 x ≠ 0 s(x)= \begin{cases} 0 &x=0\\ 1 & x\neq0 \end{cases} s(x)={01x=0x̸=0
然後定義極值函數:
m i n ( Q ( X , u ⃗ , v ⃗ , λ ) ) = m i n ( ∑ i = 1 n ∑ j = 1 p [ s ( x i j ) ∗ ( x i j − λ u i v j ) 2 ] ) s t . { ∑ i = 1 n u i 2 = 1 ∑ j = 1 p v i 2 = 1 min( Q(X,\vec u ,\vec v,\lambda))=min( \sum_{i=1}^n \sum_{j=1}^p[s(x_{ij})*(x_{ij}-\lambda u_i v_j)^2] ) \\ st. \begin{cases} \sum_{i=1}^n u_i^2 &=1\\ \\ \sum_{j=1}^p v_i^2 & =1 \end{cases} min(Q(X,u
,v
,λ))=min(i=1∑nj=1∑p[s(xij)∗(xij−λuivj)2])st.⎩⎪⎨⎪⎧∑i=1nui2∑j=1pvi2=1=1
極小化目标函數 Q ( X , u ⃗ , v ⃗ , λ ) Q(X,\vec u ,\vec v,\lambda) Q(X,u
,v
,λ),可以得到一組解 ( λ , u ⃗ , v ⃗ ) (\lambda,\vec u,\vec v) (λ,u
,v
),具體求解算法類似于 1.2節中介紹的SVD
分解的疊代算法,對于某個缺失值 x i j x_{ij} xij,隻需要計算 λ u i v j \lambda u_iv_j λuivj,就估計得到某人對某種商品的潛在需求,進而進行推薦。
2.3 代替因子分析和對應分析
未加指明時預設SVD分解中V的對角線元素是是左上到右下,按照元素平方從大到小排列的。
假定大家熟悉因子分析(不懂的自行補腦_),因子分析分為R型因子分析和Q型因子分析,對于同一組樣本,将兩種分析結合起來,就是對應分析。
2.3.1 簡單介紹因子分析模型
因子分析主要思想是将觀測名額或者觀測樣本視為某些隐藏的因素疊加而來的。比如樣本矩陣X的p個名額,分别記為 x ⃗ i , i = 1.. p \vec x_i,i=1..p x
i,i=1..p, x ⃗ i \vec x_i x
i是X的第i列。假設有r個隐藏因子記為 F j , j = 1.. r F_j,j=1..r Fj,j=1..r,建立如下模型:
x ⃗ i = u i + a i k ∗ F k + e i , ∀ i = 1... p , k = 1... r , e i ∼ N ( 0 , ε ) , ∑ i = 1 p a i k 2 = 1 \vec x_i=u_i+a_{ik}*F_k+e_i,\forall i=1...p,k=1...r,e_i \sim N(0,\varepsilon),\sum_{i=1}^pa_{ik}^2=1 x
i=ui+aik∗Fk+ei,∀i=1...p,k=1...r,ei∼N(0,ε),i=1∑paik2=1
e i e_i ei噪音, a i j a_{ij} aij是名額 x ⃗ i \vec x_{i} x
i在隐藏因子 F j F_j Fj上的得分,并且要求 c o v ( F i , F j ) = 0 , ∀ i ≠ j cov(F_i,F_j)=0, \forall i\neq j cov(Fi,Fj)=0,∀i̸=j。
表述成矩陣形式,就有:
X = ( x ⃗ 1 , x ⃗ 2 , . . . , x ⃗ p ) = U + ( F 1 , F 2 , . . . , F r ) ∗ A + e = U + F ∗ A + e , A ∈ R r ∗ p X=(\vec x_1,\vec x_2,...,\vec x_p)=U+(F_1,F_2,...,F_r)*A+\mathbb e=U+\mathbb F*A+\mathbb e,A\in R^{r*p} X=(x
1,x
2,...,x
p)=U+(F1,F2,...,Fr)∗A+e=U+F∗A+e,A∈Rr∗p
2.3.2 對比因子分析模型與SVD模型
對照這個模型,重新審視SVD,假定X的值為k,選擇 r r r滿足 r < k r<k r<k。記 u ⃗ i , i = 1... k \vec u_i,i=1...k u
i,i=1...k, u ⃗ i \vec u_i u
i是U的第i列,對應 v ⃗ i , i = 1... k \vec v_i,i=1...k v
i,i=1...k, v ⃗ i t \vec v_i^t v
it是V的第i列,$ \sigma_{i},i=1…k , , , \sigma_i 是 是 是 \Sigma 的 第 i 個 對 角 線 元 素 。 重 新 用 普 分 解 的 方 式 可 以 将 X 表 示 為 的第i個對角線元素。重新用普分解的方式可以将X表示為 的第i個對角線元素。重新用普分解的方式可以将X表示為 X = ∑ i = 1 k σ i u ⃗ i ∗ v ⃗ i t = ∑ i = 1 r σ i u ⃗ i ∗ v ⃗ i t + ∑ i = r + 1 k σ i u ⃗ i ∗ v ⃗ i t X=\sum_{i=1}^k \sigma_i\vec u_i*\vec v_i^t=\sum_{i=1}^r \sigma_i\vec u_i*\vec v_i^t+\sum_{i=r+1}^k \sigma_i\vec u_i*\vec v_i^t X=∑i=1kσiu
i∗v
it=∑i=1rσiu
i∗v
it+∑i=r+1kσiu
i∗v
it$
設 e = ∑ i = r + 1 k σ i u ⃗ i ∗ v ⃗ i t \mathbb e=\sum_{i=r+1}^k \sigma_i\vec u_i*\vec v_i^t e=∑i=r+1kσiu
i∗v
it,
對照R型因子分析,記 F = ( σ 1 u ⃗ 1 , σ 2 u ⃗ 2 , . . . , σ r u ⃗ r ) = S ′ V ′ \mathbb F=(\sigma_1\vec u_1,\sigma_2\vec u_2,...,\sigma_r\vec u_r)=S'V' F=(σ1u
1,σ2u
2,...,σru
r)=S′V′;
模型比較總結如下:記 A = V ^ t = ( v ⃗ 1 , v ⃗ 2 , . . . v ⃗ r ) t A=\hat V^t=(\vec v_1,\vec v_2,...\vec v_r)^t A=V^t=(v
1,v
2,...v
r)t,有 F ∗ A = ∑ i = 1 r σ i u ⃗ i ∗ v ⃗ i t \mathbb F*A=\sum_{i=1}^r \sigma_i\vec u_i*\vec v_i^t F∗A=∑i=1rσiu
i∗v
it,有 X = F ∗ A + e X=\mathbb F* A+\mathbb e X=F∗A+e,于是:
1、對應于因子分析中的因子得分 ∑ i = 1 p a i k 2 = 1 \sum_{i=1}^pa_{ik}^2=1 ∑i=1paik2=1,有 v ⃗ i t ∗ v ⃗ i = 1 \vec v_i^t*\vec v_i=1 v
it∗v
i=1;
2、對應于因子分析中的 c o v ( F i , F j ) = 0 cov(F_i,F_j)=0 cov(Fi,Fj)=0,有 σ i σ j u ⃗ i t ∗ u ⃗ j = { 0 i f i ≠ j u i 2 i f i = j \sigma_i \sigma_j\vec u_i^t*\vec u_j=\begin{cases}0 & if & i\neq j\\u_i^2 & if & i=j\end{cases} σiσju
it∗u
j={0ui2ififi̸=ji=j
3、對應于 e i ∼ N ( 0 , ε ) e_i \sim N(0,\varepsilon) ei∼N(0,ε),有 ( ∑ i = r + 1 k σ i u ⃗ i ) t ∗ ( ∑ i = r + 1 k σ i u ⃗ i ) = ∑ i = r + 1 k σ i 2 = ε (\sum_{i=r+1}^k \sigma_i\vec u_i)^t*(\sum_{i=r+1}^k \sigma _i\vec u_i)=\sum_{i=r+1}^k\sigma_i^2=\varepsilon (∑i=r+1kσiu
i)t∗(∑i=r+1kσiu
i)=∑i=r+1kσi2=ε
2.3.3分析總結
1. 可以看出,除了中心化的要求外,SVD分解與因子分析嚴格對應,而且SVD集R型和Q型分析兩種分析于一體,是接近于對應分析的一種分析方案。
2. 更進一步,參照普分解理論,參考因子分析,對于某個樣本的某個名額值 x i j x_{ij} xij,可以将SVD解理為: x i j = ∑ k = 1 r σ k u i k ∗ v k j + ε i j , ∀ i = 1... n , j = 1.. p x_{ij}=\sum_{k=1}^r\sigma _ku_{ik}*v_{kj}+\varepsilon_{ij},\forall i=1...n,j=1..p xij=k=1∑rσkuik∗vkj+εij,∀i=1...n,j=1..p
- σ k \sigma _k σk代表第k個隐藏因子的對名額得分的基準權重; u i k u_{ik} uik代表第i個樣本值對應第k個因子的得分; v k j v_{kj} vkj代表第j個名額對應第k個因子的得分,絕對值較小特征值代表的隐藏因子視為噪音。