2.5　正定性

正定性（positive definiteness）是凸优化问题经常出现的矩阵概念。这一小节简单了解矩阵正定性。矩阵正定性分为如下几种情况：

判断矩阵是否为正定矩阵，本册主要采用如下三种方法。

前两种方法对应代码如下：

issymmetric(A)
lambdas = eig(A)
isposd ef = all(lambdas > 0)

[~,positive_def_flag] = chol(A)
% positive_def_flag
% flag = 0, then the input matrix is symmetric positive definite

用主子式判断正定矩阵，用下例3 × 3 方阵A：

A的三个顺序主子式如下：

很容发现这三个顺序主子式均大于零，因此A为正定矩阵。

若矩阵A为负定矩阵，则A的特征值均为负值。矩阵A为半正定矩阵，则矩阵A特征值为正值或0。矩阵A为半负定矩阵，则矩阵特征值为负值或0。

下面简单说明一下矩阵特征值、特征向量和矩阵正定性关系。非零向量x为矩阵A特征向量，λ为对应特征值，则下式成立：

x为非零向量，则x^Tx > 0；若特征值λ大于0，则λx^Tx > 0，即x^TAx > 0。

A进行Cholesky分解，x^TAx写成如下形式：

R中列向量线性无关，若x为非零向量，则Rx ≠ 0，因此上式x^TAx > 0。值得注意，在一般情况，资产收益率方差-协方差矩阵都是正定矩阵，这一点在投资组合优化问题中很重要。为更直观地理解矩阵正定性，我们从几何角度来解释。

对于一个2 × 2对称矩阵A，构造如下二元函数：

在x-y-z正交空间中，当矩阵A正定性不同时，z = f（x, y）对应曲面会有不同性质，如图2.35所示。图2.35（a）所示A为正定矩阵时，z = f（x, y）曲面为凸面；该曲面又叫作椭圆抛物面（elliptic paraboloid）。当A为负定矩阵时，z = f（x, y）曲面为凹面，如图2.35（b）所示。图2.35（c）和（d）分别展示A为半正定和半负定矩阵时，z = f（x, y）曲面形状；这两个曲面分别为山谷面（valley surface）和山脊面（ridge surface）。图2.35（e）和（f）展示矩阵A为不定条件下，z = f（x, y）曲面形状，该曲面形状叫作双曲抛物面（hyperbolic paraboloid），又常被称作马鞍面（saddle surface）。本册数学部分和优化部分将会从不同角度研究这几种曲面。本节以实例讨论这几种正定性和对应曲面几何意义。

先来看一个2 × 2正定矩阵例子。矩阵A具体值如下：

容易求得A特征值分别为λ₁ = 1和λ₂ = 2，对应特征向量分别如下：

图2.36展示z = f（x, y）曲面两个不同视角视图。在该曲面边缘任意一点放置一个小球，小球都会朝着曲面最低点滚动。图2.36中点A、B和C为三个例子；曲面坡度不同，因此不同点朝下滚动时初始“加速度”大小和方向都可能会有所不同，本章后文会用梯度向量来量化该“加速度”。

再看一个2 × 2旋转正定矩阵情况。A矩阵具体如下：

经过计算得到A特征值也是λ₁ = 1和λ₂ = 2，这两个特征值对应特征向量分别如下：

z = f（x, y）曲面对应图像如图2.37。借用坐标旋转，坐标点（x, y）绕原点顺时针（clockwise）旋转θ到达（X, Y），通过下式获得：

很容易发现，（x, y）经过顺时针45°旋转到达（X, Y）点：

即

当X和Y都不为零，即[X; Y] 为非零向量时，上式恒大于0。

特殊情况下，若特征值相等，椭圆抛物面为正圆抛物面：

在图2.38曲面最小值点放置一个小球，若小球受到任何扰动，小球仍然会回落到最低点。

下面讨论一下负定矩阵情况。下式中，A为负定矩阵：

很容易求得A特征值分别为λ₁ = -2和λ₂ = -1，对应特征向量分别如下：

图2.39展示负定矩阵对应曲面；发现z = f（x, y）对应曲面为凹面。在曲面最大值处放置一个小球，小球处于不稳定平衡状态。受到轻微扰动后，小球沿着任意方向运动，都会下落。

下面来聊一聊半正定矩阵情况。矩阵A取值如下：

经过计算，矩阵A特征值分别为λ₁ = 0和λ₂ = 1；这两个特征值对应特征向量分别如下：

图2.40展示z = f（x, y）对应曲面；除了位于y轴上点以外，任意点处放置一个小球，小球都会滚动到山谷面谷底。谷底位置对应一条直线，这条直线上每一点都是函数z = f（x, y）最小值。小球在特征值为0对应特征向量方向运动，函数值没有任何变化。

下式A为半正定矩阵：

矩阵A特征值为λ₁ = 0和λ₂ = 1，对应特征向量如下：

图2.41展示旋转山谷面。同样，小球沿v₁（特征值为0对应特征向量）方向运动，函数值没有任何变化。值得指出，x和y线性相关。

下面看一个半负定矩阵情况，矩阵A具体值如下：

求得矩阵A对应特征值为λ₁ = -1和λ₂ = 0，对应特征向量如下：

图2.42展示半负定矩阵对应山脊面，发现曲面有无数个最大值。在任意最大值处放置一个小球，受到扰动后，小球会沿着曲面滚下。和山谷面一样，小球沿v₂（特征值为0对应特征向量）方向运动，函数值没有任何变化。

本节最后看一下不定矩阵情况。A为不定矩阵如下：

求得矩阵A对应特征值为λ₁ = -1和λ₂ = 1，对应特征向量如下：

图2.43展示z = f（x, y）对应曲面。当z不为零时，曲面对应等高线为双曲线；当z为零时，曲面对应等高线是两条在x-y平面内直线（图2.43中深色轨道），这两条直线即双曲线渐近线。图2.43告诉我们，曲面边缘不同位置放置小球会有完全不同结果；A点和B点处松手小球会向中心方向滚动，C点小球会朝远离中心方向滚动。

图2.43所示马鞍面中心既不是极小值点（如图2.36曲面），也不是极大值点（如图2.39曲面）；图2.43中马鞍面中心点被称作为鞍点（saddle point）。另外，沿着图2.43中深色轨道运动，小球高度没有任何变化。

图2.43中马鞍面顺时针旋转45°得到图2.44曲面。图2.44曲面对应矩阵A如下：

在z = f（x, y）为非零定值时，发现上式为反比例函数。

以下代码获得本节图像。

B4_Ch1_6.m

clc; clear all; close all
syms x y

x0 = 0; y0 = 0;
r = 2; num = 30;
[xx_fine,yy_fine] = mesh_circ(x0,y0,r,num);

A = [1, 0; 0, 1;];
% A = [1, 0; 0, 2;];
% % A = [-1, 0; 0, -2;];
% % A = [1, 0; 0, -1;];
% % A = [0, 1/2; 1/2, 0;];
% A = [1, 0; 0, -1;];
% % A = [0, 1/2; 1/2, 0;];
% % A = [1, 0; 0, 0;];
% % A = [0, 0; 0, -1;];
theta_deg = 0; % please update the theta
theta = deg2rad(theta_deg);
R = [cos(theta), -sin(theta);
      sin(theta),  cos(theta)]
inv_R = inv(R)
A_new = inv_R*A*R

f = [x, y]*A_new*[x; y];
simplify(f)
vpa(simplify(f),5)
issymmetric(A_new)
lambdas = eig(A_new)
isposdef = all(lambdas > 0)

[~,positive_def_flag] = chol(A_new)
% positive_def_flag
% flag = 0, then the input matrix is symmetric positive definite

ff_fine = double(subs(f, [x y], {xx_fine,yy_fine}));

figure(1)
c_levels =  min(ff_fine(:)):(max(ff_fine(:))-min(ff_fine(:)))/9:max(ff_fine(:));

h = mesh(xx_fine,yy_fine,ff_fine); hold on
h.EdgeColor = [1,1,1]/2;
[~,h2] = contour3(xx_fine,yy_fine,ff_fine)
h2.LevelList = c_levels;
h2.LineWidth = 1.25;
grid off; box off
hAxis = gca;
set(gca,'xtick',[])
set(gca,'ytick',[])
set(gca,'ztick',[])

xlabel('x');ylabel('y');zlabel('z');
hAxis.XRuler.FirstCrossoverValue  = 0; % X crossover with Y axis
hAxis.YRuler.FirstCrossoverValue  = 0; % Y crossover with X axis
hAxis.ZRuler.FirstCrossoverValue  = 0; % Z crossover with X axis
hAxis.ZRuler.SecondCrossoverValue = 0; % Z crossover with Y axis
hAxis.XRuler.SecondCrossoverValue = 0; % X crossover with Z axis
hAxis.YRuler.SecondCrossoverValue = 0; % Y crossover with Z axis
view(-45,60)
view(-30,90)
xlim([min(xx_fine(:))*1.2,max(xx_fine(:))*1.2])
ylim([min(yy_fine(:))*1.2,max(yy_fine(:))*1.2])
zlim([min(ff_fine(:))-1,max(ff_fine(:))+1])

function [xx,yy] = mesh_circ(x0,y0,r,num)

theta = [0:pi/num:2*pi];
r = [0:r/num:r];

[theta,r] = meshgrid(theta,r);

xx = cos(theta).*r + x0;
yy = sin(theta).*r + y0;

end