1、为什么要进行矩阵分解
个人认为,首先,当数据量很大时,将一个矩阵分解为若干个矩阵的乘积可以大大降低存储空间;其次,可以减少真正进行问题处理时的计算量,毕竟算法扫描的元//代码效果参考:http://www.ezhiqi.com/bx/art_7429.html 素越少完成任务的速度越快,这个时候矩阵的分解是对数据的一个预处理;再次,矩阵分解可以高效和有效的解决某些问题;最后,矩阵分解可以提高算法数值稳定性,关于这一点可以有进一步的说明,
借用一个上学时老师给的例子:
有方程组:
令 , ,
解方程组可得:
现在对b进行微小扰动:
,扰动项为:
此时相应的解为:
。
这个例子说明,当方程组常数项发生微小变动的时候会导致求出的结果差别相当大,而导致这种差别的并不是求解方法,而是方程组系数矩阵本身的问题,这会给我们解决问题带来很大危害,例如,我们在用计算机求解这类问题时难以避免在计算当中出现舍入误差,如果矩阵本身性质不好会直接导致所答非所问。
对常数向量b和矩阵A进行一个简单的扰动分析:
1)、扰动b,原方程组为:
(式子1),(,A非奇异)
扰动后为:
(式子2)
把式子1带入式子2得:,用2-范式来衡量这种变化得:,由于,于是得到:
而利用式子1同理可得,整理后得:
,可见b的扰动对解的影响由决定。
2)、扰动A,扰动后为:
(式子3),(,A非奇异)
稍微做一下变换:
把式子1带入后得到:
对两边同时取2-范式有:
于是有:
,整理一下就是:
//代码效果参考:http://www.ezhiqi.com/zx/art_4506.html,A的扰动对解的影响依然是由决定。
3)、对于同时扰动A和b的情况偶就不推了,最后的结果依然是,扰动对解的影响依然由决定。
定义矩阵的条件数来描述矩阵的病态程度,一般认为条件数小于100为良态,条件数在100到1000之间为中等程度的病态,条件数超过1000存在严重病态。以上面的矩阵A为例,采用2-范数表示的条件数为:,看来矩阵处于中等病态程度。
矩阵其实就是一个给定的线性变换,特征向量描述了这个线性变换的主要方向,而特征值描述了一个特征向量的长度在该线性变换下缩放的比例,有关特征值和特征向量的相关概念可查看对开篇的例子进一步观察发现,A是个对称正定矩阵,A的特征值分别为:14.93303437 和:0.06696563,两个特征值在数量级上相差很大,这意味着b发生扰动时,向量x在这两个特征向量方向上的移动距离是相差很大的——对于对应的特征向量只需要微小的移动就可到达b的新值,而对于,由于它比起太小了,因此需要x做大幅度移动才能到达b的新值,于是悲剧就发生了……………..。
关于矩阵可以有以下各种分解方式,①矩阵的三角分解(Cholesky分解、LU分解等),②矩阵的正交三角分解(QR分解等),③矩阵的满秩分解,④矩阵的奇异值分解(SVD)(关于SVD可以查看高人LeftNotEasy的一文以及他提供的参考资料)。
再看矩阵A,它是个对称正定矩阵,对这种矩阵都可以进行Cholesky分解,也就是将矩阵A分解为:,其中为一个下三角矩阵,具体操作随后讨论,回头看方程组,它就变成了:
,
将它看成两个方程组:
和,其中:,特征值为:2.2360680和:0.4472136
此时采用2-范数表示的条件数为,显然上面这两个方程组也都是良态的且只需要存储矩阵的下三角部分即可,矩阵分解的优点可见一斑。
2、实正定Hermit矩阵的完全Cholesky分解
设矩阵A有如下形式:
借用中的推导:
令,在第i次迭代时有:
,其中为i-1维单位矩阵
定义矩阵:
于是要满足,就有:
这样一直迭代下去,直到,也就是有:,最后得到分解后的下三角矩阵的元素:
j" src="">
当然也可以是上三角矩阵,此时
所以的元素是:
i" src="">
一种比较直白的分解算法可以是这样的:
设有实对称矩阵:
以及向量和保存分解结果的矩阵,算法伪码描述如下:
1: for i:=1 to n do //逐行计算L的值,w向量初始状态为实hermit矩阵A的上三角部分的一行
2: w(i .. n) := a(i,i .. n);
3: for k:=1 to i-1 do
4: temp := L(k,i);
5: if(temp != 0) then
6: w(i .. n) := w(i .. n) - temp*L(k,i .. n); //计算L(i,j)公式的分子部分
7: endif;
8: endfor;
9: w(i) := sqrt(w(i)); //L矩阵对角线元素计算
10: w(i+1 .. n) := w(i+1 .. n) / w(i); //L矩阵非对角线元素计算
11:
12: L(i,i .. n) := w(i .. n); //更新L矩阵
13: w(i .. n) //清空w向量
14: endfor;
一个简单的实现如下:
1: //串行完全Chlolesky分解,时间复杂度(O(n^3))
2: double **complete_cholesky_decompose(double **A)
3: {
4: if(NULL==A)
5: return NULL;
6:
7: double **L=malloc_matrix();
8: clear_matrix(L);
9:
10: int i,j,k,m;
11:
12: double *w=malloc_vector();
13:
14: clear_vector(w);//清除向量
15: