儆厦霞叭1. 引言复习上一篇文章《最小二乘问题详解1线性最小二乘》中的知识对于一个线性问题模型f(x;θ)Aθ那么线性最小二乘问题可以表达为求一组待定值θ使得残差的平方和最小minθ∥Aθ?b∥2本质上是求解超定线性方程组Aθb具体的线性最小二乘解是θ?(ATA)?1ATb(1)2. 求解2.1 问题虽然线性最小二乘解已经给出但是并不意味着在实际的数值计算中就能按照式(1)来进行求解。一个典型的问题就是求逆矩阵在工程实践和数值计算中直接求解逆矩阵通常是一个性能消耗大且可能不精确的操作应该尽量避免。举例来说我们按照大学本科《线性代数》课程中的方法写程序来求解一个逆矩阵假设使用伴随矩阵法A?11det(A)?adj(A)其中det(A)是矩阵A的行列式。adj(A)是A的伴随矩阵。为了求解伴随矩阵adj(A)求代数余子式 (Cofactor)对于矩阵A中的每一个元素aij计算其代数余子式Cij。代数余子式Cij(?1)ij?MijMij是删去A的第i行和第j列后得到的子矩阵的行列式称为余子式。构造余子式矩阵将所有代数余子式Cij按照原来的位置排列形成一个新矩阵C称为余子式矩阵。转置将余子式矩阵C进行转置得到的矩阵就是伴随矩阵adj(A)。adj(A)CT代入公式将det(A)和adj(A)代入公式A?11det(A)?adj(A)即可。这里我们大概能估算使用伴随矩阵法求逆矩阵的理论复杂度是O(n!)这是一个阶乘级的增长算法效率非常低。《线性代数》中介绍的另外一种算法高斯消元法也只能达到O(n3)呈指数级增加。其实效率只是一方面的问题使用计算机求解的另外一个问题是舍入误差累积在计算机中浮点数运算存在固有的舍入误差求逆过程涉及大量的除法和减法运算这些误差会在计算过程中不断累积和传播。总而言之使用通解求解逆矩阵可能存在不精确且性能消耗大的问题。2.2 QR分解那么不使用逆矩阵怎么办呢我们需要注意的是最小二乘问题的本质是求解而不是求逆矩阵因此关键是要求解正规方程ATAθATb对矩阵A作QR分解AQ1R其中:Q1∈Rm×n列正交满足QT1Q1InR∈Rn×n是上三角矩阵如果A列满秩则R的对角元均非零可逆。那么把AQ1R代入正规方程得到(Q1R)T(Q1R)x(Q1R)Tb左边整理因为QT1Q1InRTQT1Q1RxRTRx右边为RTQT1b因此正规方程等价于RTRxRT(QT1b)若R可逆即A满秩rank(A)n则RT也可逆。左右两边左乘(RT)?1得到RxQT1b.令cQ?1b这是一个长度为n的向量于是我们得到一个简单的上三角线性系统Rxc,cQT1b这就是QR方法把正规方程化简得到的核心结果只需解上三角方程RxQT1b。以上只是对A列满秩的情况做了推导如果A列满秩那么QR分解可以表示为xR?1Q?1b如果A列不满秩R奇异需要使用列主元QR方法对RTRxRT(QT1b)进行求解或者干脆使用下面要介绍的SVD分解奇异值分解法。2.3 SVD分解另外一种求解的方法是SVD分解。对任意矩阵A存在奇异值分解AUΣVT其中:U∈Rm×m为正交列为左奇异向量V∈Rn×n为正交列为右奇异向量Σ∈Rm×n为“对角块”矩阵通常写成Σ[Σr000]其中Σrdiag(σ1,…,σr)(σ1≥σ2≥?≥σr0)rrank(A)。将SVD代入正规方程先计算A?AATA(UΣVT)T(UΣVT)VΣTUTUΣVTV(ΣTΣ)VT.注意UTUI。而ΣTΣ是n×n的对角块矩阵其非零对角元就是σ2i(i1..r)其余为零。同样的计算ATbATbVΣTUTb.于是正规方程变为V(ΣTΣ)VTxVΣTUTb.两边左乘VT因为V正交VTVI得到(ΣTΣ)(VTx)ΣT(UTb)把yVTx与cUTb代入得到更简单的对角方程(ΣTΣ)yΣTc接下来按奇异值分块展开对角方程先写出Σ相关的形状Σ[Σr000],Σ?Σ[Σ2r000]对y和c也做相应分块y[y1y2],c[c1c2]其中y1,c1∈Rr对应非零奇异值y2,c2对应奇异值为0的部分维度n?r,代入得到分块方程:[Σ2r000][y1y2][Σr000][c1c2]即等价于两组方程Σ2ry1Σrc1,00?c2(无约束自由分量)由于Σr为对角且可逆第一式等价于Σry1c1?y1Σ?1rc1.而y2对应零奇异值的分量在正规方程中不受约束——这反映了在列秩不足时普通最小二乘解不是唯一的可以在零空间方向任意加解。为得到最小范数解惯常的选择取y20。最后回到x的求解对于y有y[Σ?1rc10]将c1与cU?b关系代回y[Σ?1r000]UTb由于yVTx于是xVyV[Σ?1r000]UTb定义Σ为将非零奇异值取倒数后转置得到的伪逆矩阵对角块为Σ?1r其余为0)则xVΣUTb这就是 最小二乘的 Moore–Penrose 伪逆解若A列满秩则为唯一最小二乘解由于那么ΣΣ?1SVD求解公式退化为常见的xVΣ?1UTb若秩亏它给出 在所有最小二乘解中范数最小的那个minimum-norm solution。2.4 比较从以上论述可以看到SVD分解稳定且能处理秩亏的情况但比QR分解慢复杂度高通常O(mn2)QR分解在列满秩、条件数不是太差时更快若需要判定秩或求最小范数解SVD是首选。3. 补充在最后补充一些基础知识也是笔者很感兴趣的一点那就是为什么一个矩阵A可以进行QR分解或者SVD分解呢3.1 QR分解QR分解其实非常好理解它的本质其实就是大学本科《线性代数》课程中提到的施密特Gram–Schmidt正交化。我们先复习一下施密特正交化相关的知识。设有一组线性无关向量a1,a2,…,an∈Rm我们想把它们变成一组正交再归一化后变成标准正交的向量q1,q2,…,qn。具体步骤如下取第一个向量归一化q1a1|a1|对第 2 个向量先减去在q1上的投影~q2a2?(qT1a2)q1然后归一化q2~q2|~q2|对第 3 个向量减去它在前两个正交向量上的投影~q3a3?(qT1a3)q1?(qT2a3)q2然后归一化q3~q3|~q3|一般地对第j个向量~qjaj?j?1∑i1(qTiaj)qi,qj~qj|~qj|这样得到的qi就是标准正交基且每个qj只用到了前j?1个。现在把矩阵A看成由列向量组成A[a1,a2,…,an]∈Rm×n.把施密特正交化写成矩阵形式我们得到一组正交向量Q1[q1,q2,…,qn]∈Rm×n,QT1Q1In.同时原向量可以写成ajj∑i1rijqi其中rijqTiaj把这些关系拼成矩阵形式AQ1R其中R(rij)是n×n上三角矩阵因为第j列只用到前j个qi。Q1的列正交所以QT1Q1I。3.2 SVD分解SVD分解其实也非常有意思同样也可以顺着《线性代数》中基础知识来进行推导。首先复习一下特征值和特征向量。对于一个方阵A∈Rn×n如果存在一个非零向量v∈Rn和一个实数λ使得Avλv那么λ称为 特征值eigenvaluev称为对应于λ的 特征向量eigenvector接下来复习一下什么叫做对角化。如果一个n×n矩阵A可以写成APDP?1其中D是一个对角矩阵只有对角线上有元素P是一个可逆矩阵我们就说A是 可对角化的。而且通常D的对角元素是A的特征值Ddiag(λ1,…,λn)P的列是对应的特征向量