2.4　投影

線性相關、Cholesky分解（Cholesky factorization）、特征值分解（eigen decomposition）、SVD分解（singular value decomposition）、PCA分析（principal component analysis），以及上一節介紹的矩陣線性變換等概念之間關系千絲萬縷。這一節通過投影（projection）一探究竟。

舉一個例子，主成分分析實際上尋找數據在主元空間內投影。圖2.20所示杯子，它是一個3D物體。如果想要在一張圖展示這只杯子，而且盡可能多地展示杯子的細節，就需要從空間多個角度觀察杯子并找到合適角度。這個過程實際上是將三維數據投影到二維平面過程。這也是一個降維過程，即從三維變成二維。圖2.21展示的是杯子六個平面上投影的結果。

叢書第三冊數學部分介紹過向量投影運算。投影運算一般分兩種：標量投影（scalar projection）和矢量投影（vector projection）。首先用余弦解釋標量投影，如圖2.22（a）所示，b為向量a在v方向標量投影。

下式同樣用余弦解釋矢量投影：

圖2.22（b）展示第二種計算投影方法。向量a向v投影得到向量投影a，而向量差a – a 垂直于v；據此構造如下等式：

上式中，v（v^Tv）^-1v^T常被稱作帽子矩陣（hat matrix）。帽子矩陣和最小二乘回歸有著緊密聯系，本書回歸部分會深入介紹。

如圖2.23（a）所示，兩個線性無關向量v₁和v₂構造一個平面H；圖2.23（b）所示為，多個線性無關向量v₁、v₂、…、v_q構造一個超平面F。向量a向H平面投影得到向量a。向量a由v₁和v₂構造。

向量差a – a垂直于H平面，因此垂直于平面內任何向量。

以上結論也適用于圖2.23（b）展示超平面情況。

下面來看一看數據點投影。如圖2.24所示，平面上一點Q（4, 6）和直線上不同位置點之間距離構成一系列線段，d表示這些線段長度。圖2.24下兩圖展示d和d²（線段長度平方值）隨位置變化。這些線段中最短那條，即d²最小，為Q和Q點在直線上投影點P之間距離，QP垂直于直線。尋找最短線段實際上就是優化過程。優化問題構建和求解將會在本書后文詳細介紹。

如圖2.25所示，平面上多點投影到同一條直線上，獲得一系列投影線段。當直線截距和斜率不斷變化時，這些投影線段之和不斷變化。可以想象，某個直線截距和斜率組合讓投影線段和最小。以上思路即主成分分析和主成分回歸核心。本書后面會展開講解這兩種重要分析方法。

以下代碼可完成圖2.24和圖2.25計算和繪圖，還繪制出空間多點向空間直線投影圖像。

B4_Ch1_4.m

clc; close all; clear all

v = [1;1/2];

center = [0,1];

t = -2:0.25:10;
x = v(1)*t + center(1);
y = v(2)*t + center(2);

X = [0,  0;
      1,  5;
      2, -2;
      4,  6;
      6,  0;
      8, -1;];

X_c = X - center;
b   = X_c*v/(v'*v);
X_h = b*v';
X_h = X_h + center;

fig_i = 1;
figure(fig_i)
fig_i = fig_i + 1;
plot(x,y); hold on
plot(X(4,1),X(4,2),'xb')
vectors = [x',y'] - [X(4,1),X(4,2)];
h = quiver(X(4,1)+0*vectors(:,1),X(4,2)+0*vectors(:,1)...
    ,vectors(:,1),vectors(:,2),'k');
h.ShowArrowHead = 'off';
h.AutoScale = 'off';

daspect([1,1,1])
box off; grid off

figure(fig_i)
fig_i = fig_i + 1;
subplot(2,1,1)
vector_length = vecnorm(vectors,2,2);
plot(x,vector_length)
ylim([0,9]); box off

subplot(2,1,2)
vector_length = vecnorm(vectors,2,2);
plot(x,vector_length.^2)
box off

figure(fig_i)
fig_i = fig_i + 1;
plot(x,y); hold on
plot(X(:,1),X(:,2),'xb')
plot(X_h(:,1),X_h(:,2),'xr')
plot([X(:,1),X_h(:,1)]',[X(:,2),X_h(:,2)]','k')

daspect([1,1,1])
box off; grid off

%% 3D, project points to a line in space

v = [1;1/2;2];

center = [0,1,2];

t = -2:1:4;
x = v(1)*t + center(1);
y = v(2)*t + center(2);
z = v(3)*t + center(3);

X = [0,  0,  0;
     1,  5,  2;
     2, -2,  5;
     4,  6, -1;
     6,  0,  3;
     8, -1,  1;];

X_c = X - center;
b   = X_c*v/(v'*v);
X_h = b*v';
X_h = X_h + center;

figure(fig_i)
fig_i = fig_i + 1;
plot3(x,y,z); hold on
plot3(X(:,1),X(:,2),X(:,3),'xb')
plot3(X_h(:,1),X_h(:,2),X_h(:,3),'xr')
plot3([X(:,1),X_h(:,1)]',[X(:,2),X_h(:,2)]',[X(:,3),X_h(:,3)]','k')

daspect([1,1,1])
box off; grid off

有了以上向量投影和點投影基礎，下面討論數據投影、特征值分解、SVD分解和Cholesky分解關系。圖2.26展示原始數據X，X有2列[x₁，x₂]，1000行，意味著X有兩個維度x₁和x₂，1000個觀察點。觀察圖2.26，發現x₁和x₂這兩個維度概率分布幾乎一致。經過處理，數據矩陣X已經列中心化。

回憶叢書第三冊數學部分矩陣旋轉內容。如圖2.27所示，數據矩陣X通過下式繞中心（0, 0）旋轉15°得到數據Y。

從另外一個角度來看，Z相當于X向v₁和v₂這兩個向量投影得到結果，即：

其中，

通過下列簡單計算，知道v₁和v₂這兩個向量正交。

且v₁和v₂這兩個向量為單位向量：

觀察V，發現如下等式成立：

X投影在任意正交坐標系中，該操作也常被稱作基底轉換（change of basis）。

圖2.28展示從向量v₁和v₂角度觀察數據分布情況。發現數據沿v₂方向要更為密集，方差更小；沿v₁方向更為松散，方差更大。由于數據已經中心化，矩陣Z第一列向量z₁方差即投影距離平方和平均數。

圖2.29展示數據X順時針旋轉30°得到數據Z。如圖2.30所示，發現數據Z沿v₂方向變得更為集中，方差進步一步減小；沿v₁方向數據變得更為松散，方差進一步增大。如圖2.31所示，當旋轉角度增大到45°時，圖2.32告訴我們Z沿v₂方向方差達到最小值，Z沿v₁方向方差達到最大值；并且，Z兩列數據z₁和z₂，相關性幾乎為0。這種思路即PCA分析核心。

假設數據X已經中心化，因此X樣本方差-協方差矩陣Σ_X通過下式獲得：

其中，n為X行數，即觀察點數量；注意，總體方差-協方差矩陣，分母為n。對Σ_X進行特征值分解得到：

其中，V列向量為特征向量，Λ主對角線元素為特征值。因為Σ_X為對稱矩陣。因此V列向量相互垂直。

同理，計算數據Z方差-協方差矩陣Σ_Z，并得到Σ_Z和Σ_X關系：

而對X進行SVD分解，得到：

其中：

X^TX通過下式獲得：

觀察Σ_X特征值分解和X的SVD分解結果，容易得到以下結論：

叢書第二冊隨機運動內容介紹過，通過Cholesky分解得到上三角矩陣，將相關性系數為0、服從正態分布多元隨機數轉化為服從一定相關性隨機數數據。對Σ_X進行Cholesky分解，獲得如下等式關系：

L為下三角矩陣。回憶叢書第三冊介紹的矩陣轉化內容，發現V對應旋轉操作，對應縮放操作。對Σ_X進行Cholesky分解，獲得下式：

其中，為R上三角矩陣。

若Y = [y₁, y₂] 為服從相關性系數為0標準正態分布（均值為0，方差為1）二元隨機數矩陣。那么下式獲得圖2.26中數據矩陣X：

下式驗證X方差-協方差矩陣為Σ_X：

數據Z通過R （先縮放，后旋轉V^T）獲得數據X。數據Y 和數據X相互轉換關系如圖2.33所示。而SVD分解實際上也是矩陣轉化，U和V矩陣都對應旋轉，而S對應縮放。對矩陣轉化不太熟悉讀者，參考叢書第三冊數學部分內容。

順著上述思路，我們可以把多元正態分布（multivariate normal distribution）收納到以投影為核心知識網絡中來。叢書第三冊第2章介紹多元正態分布概率密度函數矩陣表達式，如下：

其中，X =（x₁, x₂, …, x_q），代表服從多元正態分布隨機數據矩陣，每一維度隨機數為列向量（如圖2.12所示）；x為行向量，代表一個觀察點；μ_X 同樣為行向量，具體形式如下：

觀察多元正態分布矩陣表達式，發現如下看似熟悉的式子：

將上式中Σ_X替換為Cholesky分解式，得到下式：

發現上式在圖2.33旋轉（V）和縮放（）操作之前加了一個中心平移操作（x?μ_X）。由此，得到x和y關系：

以上操作正是叢書第三冊第2章中討論橢圓變形過程，如圖2.34所示。

若不考慮縮放步驟，即僅僅用旋轉和平移構造y：

得到下式：

上式取任意正數定值代表著一個多維空間橢球體。如當q = 2時，得到平面內橢圓表達式：

其中，c為任何大于0常數。

本書后面將詳細介紹更多有關橢圓和其他雙曲線內容。此外，若多元正態分布隨機數據矩陣采用X =（x₁, x₂, …, x_q）^T 形式，即每一維度隨機數為行向量，觀察點x為列向量。則多元正態分布概率密度函數矩陣表達式如下：

以下代碼獲得圖2.26～圖2.33，并且獲得特征值分解、SVD分解、PCA分析和Cholesky分解之間關系。

B4_Ch1_5.m

clc; clear all; close all

corr_rho = cos(pi/6);
SIGMA = 3*[1,0;0,1]*[1,corr_rho;corr_rho,1]*[1,0;0,1];
num   = 1000;

rng('default')
X = mvnrnd([0,0],SIGMA,num);
% R = chol(SIGMA)
% X = mvnrnd([0,0],[1,0;0,1],num);
% X = X*R;

X = X - mean(X);
theta = pi*1/12*3;
v1 = [cos(theta);
       sin(theta)];

v2 = [-sin(theta);
        cos(theta)];

V = [v1,v2];

figure(1)
plot(X(:,1),X(:,2),'.'); hold on
daspect([1,1,1])
xlim([-8,8]); ylim([-8,8]);
xlabel('x_1'); ylabel('x_2')
h = quiver(0,0,v1(1),v1(2));
h.AutoScaleFactor = 3;
h = quiver(0,0,v2(1),v2(2));
h.AutoScaleFactor = 3;
hAxis = gca;
hAxis.XAxisLocation = 'origin';
hAxis.YAxisLocation = 'origin';
box off

axes_x = [-8, 0;
             8, 0;]*V;
axes_y = [0,  8;
            0, -8;]*V;

Z = X*V;

figure(2)
plot(Z(:,1),Z(:,2),'.'); hold on
plot(axes_x(:,1)',axes_x(:,2)','k')
plot(axes_y(:,1)',axes_y(:,2)','k')
daspect([1,1,1])
xlim([-8,8]); ylim([-8,8]);
xlabel('y_1'); ylabel('y_2')
hAxis = gca;
hAxis.XAxisLocation = 'origin';
hAxis.YAxisLocation = 'origin';
box off
h = quiver(0,0,1,0);
h.AutoScaleFactor = 3;
h = quiver(0,0,0,1);
h.AutoScaleFactor = 3;

edges = [-8:0.4:8];

figure(3)

subplot(2,1,1)
histogram(Z(:,1),edges,'Normalization','probability')
xlim([-8,8]); ylim([0,0.35])
ylabel('Probability'); xlabel('y2')
box off; grid off

subplot(2,1,2)
histogram(Z(:,2),edges,'Normalization','probability')
xlim([-8,8]); ylim([0,0.35])
ylabel('Probability'); xlabel('y2')
box off; grid off

SIGMA_Z = cov(Z);
figure(4)
heatmapHandle = heatmap(SIGMA_Z);
caxis(heatmapHandle,[0 6]);
%% Conversions

[n,~] = size(X); % n, number of observations

SIGMA = (X.'*X)/(n-1)

cov(X)

[V_eig,LAMBDA] = eig(SIGMA)

[U,S,V_svd] = svd(X);

S([1,2],:)

S([1,2],:).^2/(n-1)

[coeff,score,latent] = pca(X);
% coeff, V
% score, Z
% latent, lambda

figure(5)
subplot(1,2,1)
plot(X(:,1),X(:,2),'.'); hold on
daspect([1,1,1])
xlim([-8,8]); ylim([-8,8]);
xlabel('x_1'); ylabel('x_2')
h = quiver(0,0,coeff(1,1),coeff(2,1));
h.AutoScaleFactor = 3;
h = quiver(0,0,coeff(1,2),coeff(2,2));
h.AutoScaleFactor = 3;
hAxis = gca;
hAxis.XAxisLocation = 'origin';
hAxis.YAxisLocation = 'origin';
box off

subplot(1,2,2)
plot(score(:,1),score(:,2),'.'); hold on
daspect([1,1,1])
xlim([-8,8]); ylim([-8,8]);
xlabel('z_1'); ylabel('z_2')
hAxis = gca;
hAxis.XAxisLocation = 'origin';
hAxis.YAxisLocation = 'origin';
box off

R = chol(SIGMA)
% X = mvnrnd([0,0],[1,0;0,1],num);
% X = X*R;
Z = X*R^(-1);
cov(Z)

figure(5)
subplot(1,2,1)
plot(X(:,1),X(:,2),'.'); hold on
daspect([1,1,1])
xlim([-8,8]); ylim([-8,8]);
xlabel('x_1'); ylabel('x_2')
h = quiver(0,0,v1(1),v1(2));
h.AutoScaleFactor = 3;
h = quiver(0,0,v2(1),v2(2));
h.AutoScaleFactor = 3;
hAxis = gca;
hAxis.XAxisLocation = 'origin';
hAxis. YAxisLocation = 'origin';
box off

subplot(1,2,2)
plot(Z(:,1),Z(:,2),'.'); hold on
daspect([1,1,1])
xlim([-8,8]); ylim([-8,8]);
xlabel('z_1'); ylabel('z_2')
hAxis = gca;
hAxis.XAxisLocation = 'origin';
hAxis.YAxisLocation = 'origin';
box off