2.4 投影_MATLAB金融风险管理师FRM（高阶实战）-QQ阅读男生都市网

上QQ阅读APP看书，第一时间看更新

2.4　投影

线性相关、Cholesky分解（Cholesky factorization）、特征值分解（eigen decomposition）、SVD分解（singular value decomposition）、PCA分析（principal component analysis），以及上一节介绍的矩阵线性变换等概念之间关系千丝万缕。这一节通过投影（projection）一探究竟。

举一个例子，主成分分析实际上寻找数据在主元空间内投影。图2.20所示杯子，它是一个3D物体。如果想要在一张图展示这只杯子，而且尽可能多地展示杯子的细节，就需要从空间多个角度观察杯子并找到合适角度。这个过程实际上是将三维数据投影到二维平面过程。这也是一个降维过程，即从三维变成二维。图2.21展示的是杯子六个平面上投影的结果。

图2.20　咖啡杯六个投影方向

图2.21　咖啡杯在六个方向的投影图像

丛书第三册数学部分介绍过向量投影运算。投影运算一般分两种：标量投影（scalar projection）和矢量投影（vector projection）。首先用余弦解释标量投影，如图2.22（a）所示，b为向量a在v方向标量投影。

下式同样用余弦解释矢量投影：

图2.22　向量投影和标量投影

图2.22（b）展示第二种计算投影方法。向量a向v投影得到向量投影â，而向量差a – â 垂直于v；据此构造如下等式：

上式中，v（v^Tv）^-1v^T常被称作帽子矩阵（hat matrix）。帽子矩阵和最小二乘回归有着紧密联系，本书回归部分会深入介绍。

图2.23　向量向平面和超平面投影

如图2.23（a）所示，两个线性无关向量v₁和v₂构造一个平面H；图2.23（b）所示为，多个线性无关向量v₁、v₂、…、v_q构造一个超平面F。向量a向H平面投影得到向量â。向量â由v₁和v₂构造。

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

以上结论也适用于图2.23（b）展示超平面情况。

下面来看一看数据点投影。如图2.24所示，平面上一点Q（4, 6）和直线上不同位置点之间距离构成一系列线段，d表示这些线段长度。图2.24下两图展示d和d²（线段长度平方值）随位置变化。这些线段中最短那条，即d²最小，为Q和Q点在直线上投影点P之间距离，QP垂直于直线。寻找最短线段实际上就是优化过程。优化问题构建和求解将会在本书后文详细介绍。

图2.24　平面上一点向直线投影

如图2.25所示，平面上多点投影到同一条直线上，获得一系列投影线段。当直线截距和斜率不断变化时，这些投影线段之和不断变化。可以想象，某个直线截距和斜率组合让投影线段和最小。以上思路即主成分分析和主成分回归核心。本书后面会展开讲解这两种重要分析方法。

图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.26　原始数据X统计学特点

回忆丛书第三册数学部分矩阵旋转内容。如图2.27所示，数据矩阵X通过下式绕中心（0, 0）旋转15°得到数据Y。

从另外一个角度来看，Z相当于X向v₁和v₂这两个向量投影得到结果，即：

其中，

图2.27　数据X顺时针旋转15°得到数据Z

通过下列简单计算，知道v₁和v₂这两个向量正交。

且v₁和v₂这两个向量为单位向量：

观察V，发现如下等式成立：

X投影在任意正交坐标系中，该操作也常被称作基底转换（change of basis）。

图2.28展示从向量v₁和v₂角度观察数据分布情况。发现数据沿v₂方向要更为密集，方差更小；沿v₁方向更为松散，方差更大。由于数据已经中心化，矩阵Z第一列向量z₁方差即投影距离平方和平均数。

图2.28　数据Z （X顺时针旋转15°）沿两个维度统计学特点

图2.29展示数据X顺时针旋转30°得到数据Z。如图2.30所示，发现数据Z沿v₂方向变得更为集中，方差进步一步减小；沿v₁方向数据变得更为松散，方差进一步增大。如图2.31所示，当旋转角度增大到45°时，图2.32告诉我们Z沿v₂方向方差达到最小值，Z沿v₁方向方差达到最大值；并且，Z两列数据z₁和z₂，相关性几乎为0。这种思路即PCA分析核心。

图2.29　数据X顺时针旋转30°得到数据Z

图2.30　数据Z，X顺时针旋转30°，沿两个维度分布

图2.31　数据X顺时针旋转45°得到数据Z

图2.32　数据Z，X顺时针旋转45°，沿两个维度分布

假设数据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对应缩放。对矩阵转化不太熟悉读者，参考丛书第三册数学部分内容。

图2.33　数据Y 和数据X相互转换

顺着上述思路，我们可以把多元正态分布（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所示。

图2.34　椭圆变形过程　（来自丛书第三册第2章）

若不考虑缩放步骤，即仅仅用旋转和平移构造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