BLV脉动阵列实现矩阵SVD分解

理论上，由于噪声的存在，自相关矩阵是正定的，而对于对称正定阵，它的奇异值分解等价于于特征值分解。

在之前看到的论文中，有用单边jacobi算法来求出矩阵特征值及特征向量，我在matlab上并没有实现，其中有条求theta的公式不知原作者是怎么得到的，在看了BLV阵列的文章后决定来用BLV阵列求特征值和特征向量。
BLV脉动阵列的方法如下，假设一n阶矩阵（n为偶数，若n为奇数可以加一行和一列“0”元素），将其分成(n/2)²个单元，每个单元是一个2x2阶的子矩阵。1.首先处理对角线上子矩阵，求出对应的theta_l和theta_r，并使用双边jacobi旋转（如果是对称阵则theta_l与theta_r相等）。


2.第二步，在上一步求出的theta_l和theta_r分别向行和列传递，非对角线上的二阶子矩阵在接受的角度theta后同样做双边jacobi旋转，使其对角化，如图1。
3.第三步，在所有子矩阵完成对角化后，交换元素位置，完成一轮操作；
交换规则如下：

4.第四步，判断非对角元素是否收敛，若收敛则停止算法，否则返回第一步。
BLV脉动阵列图：

其中，粗箭头表示传递角度theta，细箭头表示元素交换位置。
一开始，我也没有搞懂如何交换元素，举个例子吧，4x4阶矩阵交换元素的结果如下：

经过，这样的计算就可以求出矩阵的奇异值，但是矩阵的特征向量还没有求出，
在我看到的另一篇文章中（针对对称阵）采用了修正的BLV脉动阵列的特征向量求解方法。方法如下：
首先初始给出一个n阶单位矩阵，和求解特征值时一样，将单位阵分成2x2的子矩阵。修正的BLV脉动阵列图如下（6阶阵为例）：

其中，theta1，theta2，theta3为算法第一步中求出的对角线上子矩阵的角度theta。求特征向量的方法和求奇异值一样，首先对所有子矩阵旋转，右乘R(θ)；然后交换元素位置，完成一轮操作。交换规则如下：

下面为matlab计算结果：
输入矩阵A=[1 2 3 4;2 5 6 7;3 6 8 9;4 7 9 10];
matlab的SVD函数结果：
BLV阵列计算结果：

结果正确。
BLV脉动阵列方法优点：计算简单，并行度高，利于FPGA上实现且迭代次数少。