大数据计算方法 课件 第5-10章 矩阵的特征值与奇异值计算-基于矩阵分解的数据挖掘与分析

1大数据计算方法Chap5：矩阵的特征值与奇异值计算计算矩阵特征值的基本方法（复习）矩阵的奇异值分解计算特征值/奇异值的QR算法计算特征值/奇异值的Krylov子空间迭代法Outline2计算矩阵特征值的基本方法幂法与反幂法3特征值分解及其计算特征值计算的方法求解特征方程？计算方程中多项式系数的过程是数值敏感的高阶多项式方程(分线性方程)求根:准确度?计算量?例：

,比小一点，准确特征值为1

,但,算出两个根均为1算矩阵特征值的方法可用于解多项式方程幂法(poweriteration)求最大的特征值(主特征值)、及对应的特征向量条件：主特征值是唯一的(对实矩阵它一定是实的)是伴随矩阵Cn的特征多项式roots命令(另,NCM书第10.4节)4~7

10-9

特征值、特征向量的计算

5(PageRank算法)

每步迭代后做正交化矩阵的奇异值分解矩阵的奇异值、奇异向量奇异值分解及证明有关性质奇异值计算的敏感性6

7

到1维子空间的投影矩阵是?(singularvaluedecomposition)

0

矩阵奇异值/奇异向量

8

奇异值分解(SVD)

nn9特征值分解(r<n)(正交阵)线性代数领域的“瑞士军刀”思路:设

,

奇异值分解(SVD)证明(续)要证明的存在性根据补齐正交单位向量得到正交阵所以，存在性得证!且矩阵

唯一确定10U1为m

r列正交阵和唯一吗?非零奇异值对应的,基本上唯一矩阵奇异值/奇异向量

11

v1vrvnvr+1

u1ur

ur+1um

(m-r维)(n-r维)(r维)r维

()

rur

1u1

............奇异值分解(SVD)简化SVD分解若mn,则秩r

n,由于r常

不确定，得简化SVD为矩阵的2-范数12

------------

svd(A,'econ')

Matlab:[U,S,V]=svd(A);

s=svd(A);最大

奇异值

奇异值分解(SVD)

13注:不总成立

(F-范数下的“右逆”)

RogerPenrose,2020Nobel

PrizeforPhysics

奇异值分解(SVD)

14

+…对角阵的

奇异值是什么?

奇异值分解(SVD)eigshow演示(NCM包)2阶方阵的特征向量/奇异向量找正交的两个向量{x,y}，使得

Ax

Ay

(两对向量分别是右/左奇异向量)例子与Matlabs=svd(A);[U,S,V]=svd(A)一般，奇异值与特征值不同15思考:

得到的Ax,Ay为何总沿椭圆的轴方向?>>A=gallery(3)>>[U,S,V]=svd(A)奇异值计算的敏感性敏感性分析

为正交阵，则奇异值的绝对变动量不超过扰动矩阵的2-范数对大奇异值不敏感,对小奇异值(可能是0)仍可能很敏感例子

(奇异阵,特征值全0)>>A=gallery(5)>>formatlonge>>svd(A)(设)svd(A+eps*randn(5,5).*A)反复执行：结果：16详见NCM书第10.7节2-范数~10-11类似地分析特征值分解:

特征值/奇异值计算的QR算法基本QR迭代算法上Hessenberg约化与Givens迭代位移技术与实用的QR算法对称阵奇异值的计算非对称阵17QR算法QR迭代算法“二十世纪十大算法”之一，计算中小规模矩阵所有特征值的稳定、有效的方法实Schur分解:,Q为正交阵,S为仅含1阶/2阶对角块的分块上三角阵(拟上三角阵)QR算法是得到S的一个迭代过程：

,k=0,1,2,…基本QR算法:直接对A(k)做QR分解得到Qk不保证收敛,且每个迭代步计算量大实用的QR算法:1).先将A正交

约化为上Hessenberg阵T;2).,k=0,1,2,…,采用稀疏矩阵技巧算Qk18J.G.F.Francis和V.N.Kublanovskaya各自发明A(k)A(k+1)A(0)=AAT收敛性改善,计算效率高！(2阶块对应非实特征值)QR算法

……19000~10n3/3flops

一个QR迭代步的计算对上Hessenberg矩阵T

变执行QR迭代,得到拟上三角阵S采用Givens旋转(以4阶矩阵为例)T(k+1)仍是上Hessenberg阵!QR算法20若T(k)三对角阵,则T(k+1)也是,一个QR迭代计算量~O(n)(做QR分解),k=0,1,2,…(每个QR迭代步都是一系列正交相似变换)QR算法--实对称阵

21(做QR分解),k=0,1,2,…(变换成它需4n3/3flops)是对角阵处理实对称阵的QR算法(采用Wilkinson位移)

QR算法--实对称阵22

矩阵第i行已收敛!上Hessenberg(三对角)阵的QR迭代Wilkinson位移c为较小的常数

把Householder变换、Givens旋转依次作用到单位阵上即得到特征向量矩阵QR算法--实对称阵

23

演示程序eigsvdgui

(特征值有解析解,见NCM习题10.4)导致

QR算法--实对称阵隐式Q定理设A实对称,Q和V是两个正交阵,QTAQ和VTAV都是三对角阵.若Q的第1列与V的第1列相等,q1=v1,则Q和V本质上相等,即qi=

vi,i=2,…,n隐式位移技术一个QR迭代步Givens旋转实现QR分解,

q1由参数c和s定,它们又由24(显式位移),其中Q是对做QR分解得到c-ss

c11

xxxx11

q1再做后续旋转变为三对角可将G1直接作用于T(k)

QTAQ与VTAV本质相等确定对角元等,其他绝对值等两个场合

,其第1列由决定本质相等,

但更高效

(不存)…QR算法--奇异值奇异值的计算(A的奇异值是ATA的特征值)对A作左/右正交变换不改变奇异值

m

n用QR法求BTB(对称三对角)的特征值基于Givens旋转的隐式位移QR迭代第1个旋转阵考虑

BTB

(及可能的位移)BTB的第1列部分(双对角化)使HAW为双对角阵!

…,

且相当于做了QR迭代QR算法–非对称阵

若求特征向量矩阵,~5n3

flops

25(向后误差分析)Matlab中有关命令及其他eig命令-主要对稠密矩阵d=eig(A);返回所有特征值，使用QR算法[V,D]=eig(A);返回特征值和特征向量矩阵(可能奇异阵)

(内部处理可能会对矩阵元素进行比例调整，提高数值稳定性)可用于中小规模的矩阵eigs命令针对较大型(稀疏)的矩阵(Lanczos/Arnoldi算法)d=eigs(A,k);返回最大的k个特征值(k缺省为6)[V,D]=eigs(A,k);返回最大的k个特征值与相应特征向量hess命令:正交相似变换生成上Hessenberg阵schur命令:矩阵的Schur分解27LAPACK程序包T=schur(A);[U,T]=schur(A);相当于eig没做完还可算最小,或最接近某

值的特征值区分A是

实的或复的Krylov子空间迭代法Arnoldi过程Lanczos过程Lanczos算法非对称矩阵与SVD计算

Krylov子空间方法29比幂法收敛快(迭代步少)、能求更多特征值稀疏矩阵变得稠密!Householder变换Givens旋转矩阵向量乘可以!

Arnoldi过程30...注意：其中

可任意!...nxnnxn

wj=Avjfori=1toj

hij=viTwj

wj

=wj

-hijviend

Arnoldi-MGSv1怎么取值,得到结果与Householder约化方法一致?(隐式Q定理)

若=0,则Hk的特征值是A的特征值！

Lanczos过程--实对称阵31易证明:特征值

是T的k阶顺序主子阵，它的特征值与T的有何关系？

(差)称为Lanczos过程

Lanczos过程--实对称阵32随k增大逼近;一般对

某个k<<n,已足够准!k增大j较大时

收敛慢

Lanczos算法--实对称阵33第j个特征值的残差

迭代次数太多会使正交性丧失,因此要设上限

(Ritz值/向量)

若A实对称阵最大和最小的端特征值收敛快,中间的很慢

见testLanczos.mLanczos算法34Tk的特征值收敛到A的特征值的过程(以100阶矩阵为例)迭代步数k

Arnoldi算法/Lanczos算法35算主特征值

(n=100)

eigs命令svds命令幂法Arnoldi算法误差36大数据计算方法Chap6：优化问题与凸优化简介优化问题与凸集凸函数与凸优化问题凸优化问题举例-几何问题Outline37优化问题与凸集优化问题基础凸集38

优化问题基础39Costfunction

对称阵Lossfunction(损失函数)

优化变量最优解为最优值

根据代价函数对优化问题分类代价函数不同，对应优化问题求解难度也不同半正定二次函数是一种凸函数优化问题基础40代价函数是非负二次函数线性函数是半正定二次函数的特例如果把代价函数、以及可能存在的约束函数，都限定到

凸函数，通常能在可接受的时间内找到优化问题的最优解

凸集41

非凸集也称为凹集

(concaveset)凸包一定是凸集!

凸集42

Affinespace对称正定矩阵呢?是的凸集43

凸集凸集的和是凸集

affine凸函数与凸优化问题凸函数优化问题的一般形式可行集凸优化4445下凸函数!凸函数

凸函数46

定义域内任意直线上的函数值也形成凸函数

切面总在下方Hessianmatrix

凸函数47例子当0a1时呢?几何解释

任何向量/矩阵范数都是凸的(convexandconcave)

凸函数48

凸函数的积是凸函数吗？

凸函数49

优化问题的一般形式50实际的优化问题往往除了代价函数，还有一些约束条件要满足例子:minimax优化问题—无穷范数度量的线性拟合引入松弛变量t,问题转化为:

是这个误差的最大值!Minimizethemaximalabsoluteerrors.t.(subjectto)优化问题的一般形式51非线性优化问题的两种一般形式一个等式约束可表示为两个不等式约束minimax优化问题:s.t.s.t.s.t.(写成上面的一般形式)可行集52可行点与可行集满足所有约束条件的点称为可行点所有可行点的集合称为约束集或可行集如果一个优化问题的可行集非空，则称该优化问题可行对于带约束优化问题,即使代价函数是凸的，其求解的难度也依赖于可行集的特点s.t.(feasiblepoint)(feasibleset)好的可行集(凸集)不太好的可行集(凹集)非常不好的可行集（不连续的凹集）凸优化53

Convex

optimizations.t.

(最大化凸函数)(可行集非凸)(在凸可行集上最小化凸代价函数的问题)

凸优化54原始的minimax优化问题代价函数为虽是凸优化,代价函数不光滑,求解有难度变换后的minimax优化问题凸优化问题线性规划问题凸优化求解器MOSEK();CVX(/cvx/)例如：Minimax问题xs.t.使用高效的LP和凸优化求解算法Simplexmethod(单纯形法)Interior-pointmethod(内点法)凸优化55定理6.10:

凸优化问题的局部最优就是全局最优总结凸函数的判断方法根据定义一阶导数条件看Hessian矩阵的正定性(对足够光滑函数)降维判断的定理凸函数的非负线性组合最小化凸函数比非凸函数的最小化问题容易!术语小结56非线性优化问题,nonlinearoptimizationproblem代价函数,costfunction约束,constraint半正定矩阵,positivesemi-definitematrix半正定二次函数,positivesemi-definitequadraticfunction凸函数, convexfunction黑森矩阵, Hessianmatrix可行(约束)集,feasible(constraint)set凸集, convexset

-下水平集,

-sublevelset凸优化问题, convexoptimizationproblem线性规划, linearprogramming凸优化问题举例-几何问题几何优化问题求最大内切椭球57许多几何问题可转化为凸优化问题我们以求最大内切椭球的问题为例推导出它对应的凸优化问题几何问题58求最大内切椭球求最小外接椭球问题定义给定一个多面体求其内部的体积最大的椭球椭球的表示椭球内点的集合求最大内切椭球5901

坐标放缩、旋转、平移

d

x1x2

求最大内切椭球60

两个

特征值为

求最大内切椭球61

上确界

这就是优化变量须满足的约束条件最大内切椭球

求最大内切椭球62

正交变换相当于坐标轴旋转,面积不变

椭圆方程

高维情况也成立

求最大内切椭球63

还缺什么吗?

对称正定阵Z,存在对称正定阵

使得对称阵的特征值

0为凹函数!

求最大内切椭球64

它是一个凸优化问题!求最大内切椭球的问题可形式化为:65大数据计算方法Chap7：数据分析中的凸优化最大间隔分类与支撑向量机线性回归的过拟合与正则化压缩感知技术简介Outline66最大间隔分类与SVM二分类问题支撑向量机方法(线性SVM)更多讨论67

分类问题68

更多高维数据的分类问题Bishop,PatternRecognitionandMachineLearning,2007

支撑向量机方法69

设正、负支撑面的方程分别为:

支撑向量机方法70

代价函数为凸函数，线性约束；是凸优化问题!

支撑向量机方法71

根据右图看出，最优解为：

更多讨论72

正则化参数凸优化!

线性回归的过拟合与正则化线性回归的过拟合正则化-岭回归正则化-LASSO73过拟合问题74也反映模型复杂程度明显的过拟合！

过拟合问题75例7.4评价拟合模型的误差用函数的范数如果无法增加数据点、

模型复杂度也不能降

低，如何避免过拟合？

拟合误差随多项式次数变化的情况正则化76对例7.4的更多观察各个拟合多项式的系数让多项式系数变小，则控制了它偏离准确模型的幅度对优化问题做点修改：惩罚因子

的取值：交叉验证、gridsearch当次数K很大时，拟合系数

普遍很大！在数据数量、模型复杂程度不变情况下，可否限制拟合系数的值不能太大？惩罚项Ridgeregression,Tikhonov正则化正则化77确定最优惩罚因子

的交叉验证方法将已有的数据分的一部分做为训练集、另一部分做验证4-折交叉验证(4-foldcrossvalidation)3/4数据训出(拟合出)模型另外1/4算模型误差计算误差的均值，作为

某个

对应的准确度对一系列

值，重复上述计算，选出最好的

值这是在没有准确拟合解的情况下的一个选择最优参数的方法正则化78“岭回归”问题的求解相当于要求解线性最小二乘问题

>>At=[A;sqrt(lam)*eye(n)];>>x=At\[b;zeros(n,1)];正则化79例7.5用正则化做9次多项式拟合设置不同的

值看拟合函数误差9次多项式的模型矩阵与低次多项式拟合准确度

差不多

模型足够复杂、又不过拟合

正则化-LASSO80

或

正则化-LASSO81

压缩感知技术简介图像恢复的压缩感知技术OMP算法82图像恢复的压缩感知83图像修复问题已知部分像素点的值,如何恢复整个图像?灰度图

矩阵，彩色图3个矩阵图像压缩感知技术将数据变换到频域(离散余弦变换)利用信号的稀疏性通过求解欠定的线性方程组实现信号的重建或恢复10%的像素频域数据一般是稀疏的！

84图像恢复的压缩感知

85图像恢复的压缩感知

求解思路1：松弛为凸优化问题求解思路2：用启发式方法求近似最优解86正交匹配追踪算法

87正交匹配追踪算法

这就是正交匹配追踪算法Orthogonalmatchingpursuit(OMP)88正交匹配追踪算法算法描述

主要计算量时间复杂度89正交匹配追踪算法例7.10(求解带模型稀疏约束的曲线拟合)用OMP算法依次选出的模型是：0-范数正则化不能解决过拟合问题质量不好三次多项

式误差19%过拟合现象90大数据计算方法Chap8：非线性方程与无约束优化求解求解非线性方程的方法简单的无约束优化方法基于求导的无约束优化方法Outline91求解非线性方程的方法二分法牛顿法割线法逆二次插值法解非线性方程组的牛顿法92非线性方程基本理论

93

二分法

94

0

小结

95>>f=@(x)x-1-1/x;牛顿法原理利用曲线上点(xk,f(xk))处的切线,求xk+1计算公式程序my_Newtonx*xk+1

xkxy(xk,

f(xk))

k=0;whileabs(x-xprev)>tol*abs(x)xprev=x;x=x-f(x)/fprime(x);k=k+1;end

若对于收敛快的算法，是误差的估计若tol=eps,基本上达最准其他优点：可以算复数根、解方程组96若准确解为0?修改循环结束条件(f(x)=0时退出循环)牛顿法的问题

97对任意初值都不收敛的例子一阶导数间断点另两个缺点：2阶收敛速度!割线法

x*xk

xk-1xy

k=0;whileabs(b-a)>eps*abs(b)c=b;%b是x的近似解b=b+(b-a)/(f(a)/f(b)-1);a=c;%a是上一个近似解k=k+1;end

超线性收敛速度：

98缺点:需要两个较好的初值逆二次插值法

抛物线xk+10xyxk

xk-1

xk-2

也可用polyfit,polyval算99

逆二次插值法插值函数P满足:则是

它与x轴交点程序优缺点收敛速度快f(xk-2),f(xk-1),f(xk)虽不等,若值很接近,

得到的xk+1将远远偏离求解区间割线法也有这种不稳定性(斜率~0)k=0;whileabs(c-b)>eps*abs(c)x=polyinterp([f(a),f(b),f(c)],[a,b,c],0)a=b;b=c;c=x;k=k+1;endIQI算法(Inversequadraticinterpolation)100

xf(x)

...

通用求根算法zeroin算法的特点每步迭代都使有根区间缩小抛弃逆二次插值/割线法得到的”不满意”解按”逆二次插值,割线法,二分法”的优先顺序生成下一步解，保证较快的收敛速度算法稳定、通用性强，是Matlab命令fzero的基础NCM中的演示程序fzeroguiNCM中的fzerotx.m101>>fzerogui>>bessj0=@(x)besselj(0,x);>>fzerogui(bessj0,[03.83])缺省f(x)=x^3-2x-5fzero命令

>>bessj0=@(x)besselj(0,x);>>ezplot(bessj0,[0,10*pi]);forn=1:10z(n)=fzero(bessj0,[(n-1)n]*pi);endholdon;plot(z,zeros(1,10),'o')102optimset设置结构(求单变量连续函数的零点)>>opt=optimset('tolX',1e-2);非线性方程组

103

非线性方程组例:用牛顿法求解104

解方程

解方程

该方程的Jacobi矩阵

简单的无约束优化方法黄金分割搜索下山单纯形法非线性最小二乘问题105无约束优化

106

不涉及导数

(本节)涉及导数

(下一节)例如：黄金分割搜索

*abuv

黄金分割搜索！

107少算一次f(x)黄金分割搜索进一步改进黄金分割搜索很稳定,收敛速度仍慢类似割线法/IQI法,用抛物线近似被优化函数的局部(需3个点做插值,如其函数值”高-低-高”时)何时选用抛物线的解?

在区间内,相邻解间距缩小,等fminbndMatlab命令fminbnd,求单变量函数区间上最小值[x,fval,exitflag,output]=fminbnd(fun,x1,x2,opt)x1<x2,为区间两端点;fval为目标函数(极小)值exitflag退出标记,output包含算法、函数求值次数等用optimset设置opt中选项的值(同fzero)可再与黄金分割搜索结合！abuv108*(类似于zeroin算法,见fminQI程序)fminQI程序程序停止的判据黄金分割搜索与抛物线法whileabs(x-xm)>tolr=(x-w)*(fx-fv);q=(x-v)*(fx-fw);p=(x-v)*q-(x-w)*r;q=2.0*(q-r);ifq>0.0,p=-p;endq=abs(q);r=e;e=d;%Istheparabolaacceptable?para=((abs(p)<abs(0.5*q*r))&(p>q*(a-x))&(p<q*(b-x)));ifparad=p/q;end含极小值

区间的中点x+p/q在区间[a,b]内,且相邻解间距缩小根据x,w,v这三点数据造抛物线p,或q取相反数109fminbnd的简化版本u=fminQI(F,a,b,tol,varargin)近似极小点详见课本8.2.1小节的fminQI程序最优化与fminbndfminbnd/fminQI的例子求-humps(x)的极小值搜索过程显示(习题8.5)110

>>F=@(x)-humps(x);>>fminQI(F,-1,2,1.e-4)fminbnd缺省的tolX是1e-4但结果与fminQI的不一样设opt=optimset('tolX',1e-6)试试?

下山单纯形法主要思路黄金分割搜索和fminbnd只能用于单变量问题在多维空间,通过计算单纯形顶点处函数值，然后再对它做调整(反射、扩展、收缩等)，不断重复操作，找到最小值点也叫Nelder–Mead法，是Matlab命令fminsearch的基础n维空间的单纯形n+1个非退化的点构成的形体111(2维空间的例子)3个不共线的点构成的三角形是一个单纯形这3个点共线,不构成单纯形

下山单纯形法112二维空间的例子最差点新点单位阵的第i列最好点扩展操作(expansion)向代价函数最小

顶点的方向扩展收缩操作(contraction)如果单纯形的形状

比较“狭长”两个顶点向代价函

数变小的方向收缩收缩操作2(shrink)扩展操作的逆缺点:对高维的问题收敛速度慢、可能不收敛、局部极小优点:不需要算导数(代价函数不光滑...)下山单纯形法113最好点fminserach命令求多元实函数的最小值，使用下山单纯形法[x,fval]=fminsearch(fun,x0,options)可分离的非线性最小二乘问题拟合公式中含非线性参数与线性参数例:求解最小二乘问题若已知非线性参数(

1/2)的值,则用\可求线性参数值(

1/2)最小化的是残差的2-范数:它是非线性参数(

1/2)的函数非线性最优化+线性最小二乘fminserach命令及其应用

114后面加给fun的参数带入fminsearch中的fun比全看成非线性参数直接用fminsearch求解更有效/可靠例子数据:t=(0:.1:2)'y=[5.89553.56392.51731.9790...,0.13700.22110.17040.2636]';

1,2的初值为:[36]'求解过程与演示线性最小二乘的解定义为函数:

res=expfitfun(lambda,t,y)非线性最小化:fminsearch(@expfitfun,lambda0,[],t,y)程序细节与动态演示程序

expfitdemofminserach命令及其应用(放射性物质衰变的观测)最小化过程收敛时对应的拟合曲线将观测数据设为函数参数>>editexpfitfun.m;>>editexpfitdemo.m;115基于求导的无约束优化方法梯度下降法牛顿法共轭梯度法116无约束优化方法分类黄金分割搜索及其变种算法下山单纯形法(downhillsimplexmethod)梯度下降法(gradientmethod)牛顿法(Newtonmethod)共轭梯度法(conjugategradientmethod)梯度下降法若代价函数光滑,其导数

可用于搜索最优解117

不涉及导数涉及导数

(本节)另一个思路:解方程无约束优化

118

梯度为0解单变量

无约束优化无约束优化

119

非线性方程组,用牛顿法解

=0

代价函数减少！

对称正定近似计算Hessian矩阵或它的逆无约束优化

120

无约束优化

121

Fletcher-Reeves公式

Polak-Ribiere公式存在一些经验性的尝试

Polak-Ribiere公式

如果收敛停滞,每N次迭代后可重置搜索方向为负梯度方向理论上严格下降,且收敛无约束优化涉及导数的方法小结对一些实际的高维空间优化,梯度下降法与共轭梯度法往往是比较好的选择122

fminsearch适合仅含几个或十几个变量的小规模优化问题Matlab优化工具箱(Toolbox)无约束极小化fminunc最小二乘问题也属于最优化问题特别针对最小二乘问题的函数:lsqnonlin,lsqcurvefit带约束极小化fminconMatlab中的多变量最优化解法123(包含内点法等多种算法)(包含拟牛顿法等多种算法)124大数据计算方法Chap9：约束优化解法与对偶原理等式约束优化问题不等式约束优化问题对偶方法与KKT条件Outline125求解等式约束优化问题线性等式约束子空间约减方法拉格朗日乘子法(重点：线性等式约束的两种解法)126带约束的非线性优化

127s.t.

s.t.子空间约减法

128s.t.pxn

p列n-p列

非奇异上三角方阵

s.t.

nx(n-p)

p列p列

子空间约减法129s.t.pxn做LU分解的开销小于对A做QRCP分解即使A为稀疏阵,变量代换后优化变量前的矩阵F都是很稠密的,将导致后续求解无约束优化问题时开销大一般的等式约束问题

拉格朗日乘子法130s.t.

曲面的法向为

即垂直于这p个切面的交集拉格朗日乘子法可处理一般的等式约束问题等价于例子:131s.t.s.t.最优解

垂直于g(x)等值面,即沿其梯度方向可用牛顿法解方程组未必好解；未必是凸优化极值点x*满足:

拉格朗日乘子法

132s.t.

,其中pxnn个方程求解方法：

拉格朗日乘子-牛顿法

拉格朗日乘子法处理线性等式约束若代价函数对一般的非线性代价函数,将方程整体看成用牛顿法求解迭代公式:133s.t.pxn比子空间约减法易保持稀疏性任意初始,k1后,每步解此方程还可以使用阻尼牛顿法等技巧处理线性等式约束例子:用牛顿法求解设拉格朗日乘子法134s.t.pxn可行解

s.t.若不是可行解若已是可行解求解不等式约束优化问题线性等式约束+不等式约束内点法—对数障碍法半正定规划135主要考虑凸优化假设等式约束均为线性等式约束注意:一般情况下，可行集未必是凸集处理不等式约束定义一个指示函数(indicator)注意:

在计算过程中也要保证新的优化问题只含线性等式约束,但函数I()不可导,无法直接使用拉格朗日乘子-牛顿法带不等式约束的优化问题136s.t.s.t.要求,否则代价函数无穷大u0I(u)凸函数处理不等式约束用带参数的光滑函数近似函数I()相继求解一系列问题,其中每个都

可以用牛顿法可靠地求解对数障碍法(logarithmicbarrier)带不等式约束的优化问题137s.t.s.t.homotopymethod参数t

另一种障碍函数当t0时非常近似I(u)障碍函数对数障碍函数消去不等式约束,且代价函数光滑注意:计算中仍要保证hj(x)<0例:新的代价函数是凸的吗?内点法138s.t.s.t.x*(t)处等值面与超平面相切从可行集内部的点开始也是凸函数!s.t.若是凸函数,半正定

对数障碍法1.选择初始t值和初始解x(0)2.求解线性等式约束的非线性优化

得到最优解x*

(拉格朗日乘子-牛顿法)3.令x(0)=x*,t=t(取值10~20)4.重复上述第2、第3步,直到x*趋于收敛如何选初始的x(0)

?必须是可行的，即预处理步骤:求解整体流程:对数障碍法解预处理问题得到初始解,再解原问题内点法139s.t.s.t.凸优化凸函数注意：得到的是近似解，总是不满足

约束中的=s.t.若最优解s*>0,原问题不可行;否则,得到可行的初始解用对数障碍法解它初始解容易设置内点法140s.t.凸函数特殊的不等式约束例:最大内切椭球问题这种广义不等式表示矩阵半正定,该问题称为半正定规划半正定规划一般形式:处理半正定约束条件的方法一定是凸函数半正定规划141(Semidefiniteprogramming)s.t.有些约束条件不是此形式

可看成是广义的不等式s.t.对称矩阵(其他常见约束)上面的属于这一般形式吗?

使用相应的对数障碍函数:线性函数与凸函数的复合可行集凸吗?

用对数障碍法求解要保证正定怎么选初始x(0)是问题在预处理步骤,求解解预处理问题时,

可任意取初始x(0),总可以找到足够大s(0)使得满足约束小结半正定规划s.t.s.t.若最优解s*>0,原问题不可行,否则得到可行的x使用对数障碍法线性等式约束二次优化线性等式约束优化

含不等式约束的优化

(包括半正定规划)解一次线性方程组拉格朗日乘子-牛顿法内点法(对数障碍法)Ipopt,/coin-or/Ipopt142

对偶方法与KKT条件对偶函数对偶问题强对偶KKT条件143带约束非线性优化问题拉格朗日函数:它是关于x的非线性函数,但关于u和v是线性函数拉格朗日对偶(dual)函数拉格朗日对偶函数144s.t.拉格朗日乘子

一定是个凹函数!(下确界)若x*是优化问题的解且v

0,则定理9.2：d(u,v)是f(x*)的下界例9.8:拉格朗日对偶函数145s.t.拉格朗日对偶函数:(当v

0时)s.t.

例9.7:可行集为区间[-0.46,0.46]f(x)周围的虚线表示拉格朗

日函数

在v=0.1,0.2,…,1时的曲线拉格朗日对偶函数当0v1时,d(v)的曲线?且d(v)是凹函数拉格朗日对偶函数146s.t.,当v>0时f(x),h(x)如右图对偶问题:例9.9:对偶问题一定是凸优化问题!对偶问题147s.t.原问题s.t.v是对应的拉格朗日系数s.t.对偶问题:s.t.其解u*,v*满足：称为弱对偶条件为最优对偶间隔(optimaldualitygap)

s.t.即弱对偶普遍成立通过解对偶问题,可得原问题最优值的下界强对偶:并非总成立;若原问题为凸优化,经常成立对凸优化问题若可行集有内点,即可行集中

使得,

则强对偶成立。精细的Slater条件:若某些hj是线性函数,则可以去掉条件中它们对应的<0要求强对偶148s.t.对偶问题s.t.s.t.,均为凸函数这个条件称为Slater条件这些都是判断强对偶的充分条件!即最优对偶间隔为0定理9.4线性规划问题都满足Slater条件！

强对偶149对偶问题可表示为:凸的半正定规划!非凸优化、但满足强对偶对偶函数由于互补松弛(complementaryslackness)若强对偶成立,x*为原问题最优解,u*,v*为对偶问题最优解等号成立,则,j=1,2,…,qKKT条件(Karush-Kuhn-Tucker)若强对偶成立,且相关函数可微

是的最小值点KKT条件150s.t.注意:(互补松弛条件)若则若则激活约束未激活约束(驻点条件)

KKT条件151s.t.

KKT条件与强对偶都成立问题(a):小结152对于强对偶成立的凸优化问题,KKT条件也是充分的优化问题的解法总结内点法+牛顿法(主要对凸优化)如果强对偶条件成立(如针对凸优化的Slater条件,…),求解对偶问题有时更简单（它一定是凸优化）可基于KKT条件求解（一组方程+不等式）无约束优化线性等式约束优化不等式约束优化(线性等式)、半正定规划(也适合非线性等式约束)黄金分割搜索及其变种下山单纯形法梯度下降法牛顿法共轭梯度法子空间约减法拉格朗日乘子-牛顿法KKT条件满足并不意味着存在强对偶153大数据计算方法Chap10：基于矩阵分解的数据挖掘与分析PCA与降秩最小二乘随机化矩阵低秩分解非负矩阵分解Outline154PCA与降秩最小二乘PCA不满秩最小二乘问题主分量回归155主成分分析(PCA)数据矩阵特征间有相关性(数据有误差),可用少量独立特征(主成分)刻画找一组特征的主成分z1…,zk,使n个特征近似为它们的组合156第1次观测数据第m次观测数据第i

次观测数据特征1特征nA秩

k主成分

即A的前几个左奇异向量怎么求主成分?

normalizedfirstprincipalcomponent低秩近似(分解)主成分分析(PCA)一组高度与重量的数据两列数据是非常有关联的近似表示高度与重量有一个潜在的主要成分>>[U,S,V]=svd(A,0);>>sigma=diag(S)sigma= 156.4358 8.7658>>E1=sigma(1)*U(:,1)*V(:,1)'第一主成分是：157V(:,1)'=[0.9468,0.3219]主成分方向principalcomponentdirections主成分分析(PCA)

158

(居中化)主成分分析(PCA)非负矩阵的第一对奇异向量矩阵A的各行数据到其距离平方和

最小的过原点直线,一定如图穿过

这堆数据点根据之前的结论,该该直线方向沿

右奇异向量v1所在方向若A为非负矩阵,则向量v1中元素全非负,或全非正若取v1为非负向量,则,即左奇异向量u1也非负对非负矩阵,其第一个左/右奇异向量可同时为非负向量!PCA用途降维(分类、聚类的预处理,可视化),嵌入表达,LSI,…1590101010(第一主方向)若A还对称呢?主特征值一定是正数通常要先对数据矩阵做“列居中”AA-

(sumofentriesineachcolumniszero)

降秩最小二乘

(线性回归的最小范数解)Ax

b>>x=pinv(A)*b;最小化要求是范数最小的最优解解求最小范数解比基本解计算量大些2-范数相同最小范数解Rank(A)=r,列不满秩(SVD分解)(QRCP分解)160降秩最小二乘(ReducedRankLS)Ax

b由于数据误差,精确求解往往不好很多时候，x代表模型，或是中间结果分类问题:测试数据b,关心或的值回归分析:161特征1

特征nAb因变量解得到，代表模型若A的各列近似相关,这个问题很敏感新观测怎样让预测结果准确/合理?应将A近似为一个低秩模型再做拟合！overfitting第1次观测数据第m次观测数据第i

次观测数据降秩最小二乘(ReducedRankLS)

162计算很准确>>A0=[101111];>>A=[A0,A0*[1;0.5]+1e-7*randn(3,1)];>>b=A*[1;1;1]+1e-4*randn(3,1);>>x=A\b

与扰动差不多这个结果好得多!

抑制过拟合降秩最小二乘(ReducedRankLS)Ax

b在某个子空间Zk中搜索最优x,该子空间基为{z1,z2,…,zk}求解主成分回归如果A是低秩矩阵，则精确(最小范数)解x是{v1,v2,…,vr}的线性组合163解残差平方为即对A做秩k近似再求解!基于截断SVD的近似保留了A中数据的主成分!取前k个右奇异向量形成Zk?主成分回归主成分回归(PrincipalComponentRegression)选择适当的k,解例10.4:文档-关键词矩阵查询向量,与文档库的匹配度q1:“排名”,“网页”,“万维网”q2:“中国队”,“国际足联”164若解原始LLS问题,即k=5,两个残差区别不大但显然q1与文档库内容更相符若考虑k=1~3,两个残差相差很大；例如k=2,3时…解q1q2应用:文档自动归类？随机化矩阵低秩分解子空间嵌入固定秩参数的随机化低秩分解自适应的随机化低秩分解少遍历与单遍历算法有关应用与发展165子空间嵌入(subspaceembedding)

166接近最优，时间短/易并行=

m

nm

k

k

n

m

nn

l标准高斯m

l或

(离散余弦变换)(随机取列)(1/-1)固定秩参数的随机化低秩分解

167=

正交化

m

nm

l

l

n

幂迭代/同步迭代

奇异值衰减更快正交化减少数值误差

~

固定秩参数的随机化低秩分解

168

m

nm

l

l

n

~

l

l

l

l

l

n近似的k-SVD这可以提高结果的近似程度

（无幂迭代）

固定秩参数的随机化低秩分解基础的随机化SVD算法例10.5169>>U=orth(randn(4000,2000));>>V=orth(randn(2000,2000));>>S=diag(1./(1:2000).^2);>>A=U*S*V';>>tic;[Us,Ss,Vs]=svds(A,100);toc;历时2.569249秒。>>tic;[Ub,Sb,Vb]=basic_rSVD(A,100,105,2);toc;历时0.092135秒。>>norm(diag(Ss)-diag(Sb))ans= 8.0063e-06奇异值为的矩阵构造奇异值按平方规律衰减的矩阵,比较随机化算法与svds快28倍计算出的奇异值的误差很小

自适应的随机化低秩分解

170

m

nm

k

k

n

我们发现残差的Frobenius范数可以增量计算!

[1]P.-G.Martinsson,S.Voronin,Arandomizedblockedalgorithmforefficientlycomputingrank-revealingfactorizationsofmatrices,SIAMJ.Sci.Comput.,2016.svds是标准做法；随机化算法randQB在准确度要求不高、要求数据遍历少的场景下有很大优势W.Yu,Y.GuandY.Li,"Efficientrandomizedalgorithmsforthefixed-precisionlow-rankmatrixapproximation,"

SIAMJournalonMatrixAnalysisandApplications,39(3):1339-1359,Aug.2018自适应的随机化低秩分解randQB_EI算法(cont’d)171

m

nm

k

k

n数学上与randQB算法等价适合于稀疏的A节省内存用量与计算时间

>>helpsvdsketch

自适应的随机化低秩分解randQB_EI的实验结果172

自适应的随机化低秩分解randQB_EI的实验结果(不用看红色的线)173w/opowerschemew/powerscheme对于稀疏矩阵(0.3%稠密度)，比randQB_b算法最多快22倍，节省79倍内存对于稠密矩阵，有超过2倍的时间与内存节省自适应的随机化低秩分解randQB_EI的自适应分解风景图片(3168x4752)10%近似准确度7X矩阵元素压缩,比SVD快12X对用TD-IDF模型得到的关键词/

作者矩阵(8300x100000)，也得

到近似最优k值,比svds快10X

更快的farPCA算法174randQB_EISVDRangeFinderimageP=2441AMinerP=1244021158018P=22229确定出的秩k

少遍历与单遍历算法

175DDR4内存：~50GB/sSDD硬盘：~3GB/sSATA硬盘:~300MB/s

少遍历与单遍历算法PerSVD算法动态位移技术，

提高近似质量用奇异值分解代

替QR做正交化用特征值分解来

加速计算如果p=0，则得

单遍历PCA算法176使用加速技巧的动态位移技术处理硬盘上大的稠密矩阵，达到很好的速度-准确度平衡少遍历与单遍历算法任意流单遍历算法在流数据场景，PerSVD(p=0)假设数据逐行到达，未必成立实际上，可能每次数据更新为任意位置矩阵元素(任意流)Tropp算法:用随机投影建三个草图矩阵[1]JTropp,etal.,“Streaminglow-rankmatrixapproximation...,”SIAMJ.Sci.Comput.,2019

有关应用与发展dashSVD+奇异值门限(SVT)算法，用于矩阵补全利用信息冗余/低秩SVT算法使用dashSVD算法178迭代求解过程：多次截断

奇异值分解SVT+dashSVD

SVT+rSVD-BKI1218在含两个8核IntelXeonCPU@2.10GHz的Ubuntu

服务器上

SVT+svds(算法8)8有关应用与发展randQB_EI与协同过滤结合用于评分估计User/item矩阵A,基于SVD的协同过滤(CF)评分估计:计算item之间相关性,已知数据加权求和隐式因子个数k如何自动选取?randQB_EI算法+验证集数据稀疏矩阵rSVD的加速技巧测试结果接近最优,计算量非常小179(item隐式因子矩阵)最优有关应用与发展超大规模网络的嵌入(networkembedding)NetMF算法核心挑战在于第2、第3步Tropp算法+sparse-sign随机嵌入可以实现不构建M矩阵，

将计算复杂度降为线性(对节点数n、边数m)实验结果:对网络