双边Jacobi方法求解任意矩阵的奇异值分解的Matlab代码

以下是双边Jacobi方法求解任意矩阵的奇异值分解的Matlab代码：

function [U, S, V] = bidiagonal_svd(A)
% 双边Jacobi方法求解任意矩阵的奇异值分解

[m, n] = size(A);

% 初始化
U = eye(m);
V = eye(n);

% 双边Jacobi迭代
while 1
    % 判断是否达到精度或迭代次数是否超过限制
    if max(max(abs(tril(A, -1)))) < eps || max(max(abs(triu(A, 1)))) < eps
        break;
    end
    
    for i = 1:n-1
        % 计算Givens旋转角度
        [c, s] = givens(A(i, i), A(i+1, i));
        
        % 对A进行Givens旋转
        G = [c s; -s c];
        A(i:i+1, :) = G * A(i:i+1, :);
        V(:, i:i+1) = V(:, i:i+1) * G';
    end
    
    for i = 1:m-1
        % 计算Givens旋转角度
        [c, s] = givens(A(i, i), A(i, i+1));
        
        % 对A进行Givens旋转
        G = [c s; -s c];
        A(:, i:i+1) = A(:, i:i+1) * G;
        U(:, i:i+1) = U(:, i:i+1) * G';
    end
end

% 提取奇异值
S = diag(diag(A));

end

function [c, s] = givens(a, b)
% 计算Givens旋转角度

if b == 0
    c = 1;
    s = 0;
elseif abs(b) > abs(a)
    t = -a / b;
    s = 1 / sqrt(1 + t^2);
    c = s * t;
else
    t = -b / a;
    c = 1 / sqrt(1 + t^2);
    s = c * t;
end

end

输入任意矩阵A，输出其奇异值分解的左奇异向量矩阵U、奇异值矩阵S和右奇异向量矩阵V。其中，Givens旋转用于将A变为上（或下）三角形矩阵，提取奇异值后即为S，左奇异向量矩阵U和右奇异向量矩阵V分别为做、右上（或下）三角形矩阵的单位正交矩阵。