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