基于双边Jacobi的奇异值分解的Matlab代码

以下是基于双边Jacobi的奇异值分解的Matlab代码，其中A为输入的矩阵，U、S、V为输出的分解结果：

function [U, S, V] = bidiag(A)
% 双边Jacobi的奇异值分解
[m, n] = size(A);
U = eye(m);
V = eye(n);
B = A;
for i = 1:n
    for j = 1:i-1
        % 双边Jacobi旋转
        [B, Q] = jacobi(B, j, i, 'b');
        U(:, [j i]) = U(:, [j i]) * Q;
    end
    for j = i+1:n
        % 双边Jacobi旋转
        [B, Q] = jacobi(B, i, j, 'b');
        V(:, [i j]) = V(:, [i j]) * Q;
    end
end
S = diag(B);
end

function [B, Q] = jacobi(A, p, q, flag)
% 双边Jacobi旋转
n = size(A, 1);
theta = (A(q,q)-A(p,p))/(2*A(p,q));
if theta >= 0
    t = 1/(theta + sqrt(1+theta^2));
else
    t = -1/(-theta + sqrt(1+theta^2));
end
c = 1/sqrt(1+t^2);
s = c*t;
Q = eye(n);
Q(p,p) = c;
Q(q,q) = c;
Q(p,q) = s;
Q(q,p) = -s;
if flag == 'b'
    B = Q'*A*Q;
else
    B = A*Q;
end
end

代码实现了双边Jacobi的奇异值分解过程，其中包括了双边Jacobi旋转的实现。具体而言，对于输入的矩阵A，首先初始化U和V为单位矩阵，然后对于每一列i，将列i之前的所有列j（j<i）对应的双边Jacobi旋转应用于A和U，将列i之后的所有列j（j>i）对应的双边Jacobi旋转应用于A和V，最后得到分解结果U、S、V。其中，S是对角矩阵，其对角线上的元素为A的奇异值。