基于双边jacobi的svd算法Matlab代码

以下是基于双边Jacobi的SVD算法的Matlab代码：

function [U,S,V] = bidiag_svd(A)
% 双边Jacobi SVD算法
% 输入参数：
% A - 待分解的矩阵
% 输出参数：
% U,S,V - SVD分解结果

[m,n] = size(A);
B = A;
U = eye(m);
V = eye(n);

for i = 1:n
    for j = (i+1):n
        [G,R] = planerot(B(:,i),B(:,j));
        B(:,i) = G'*B(:,i);
        B(:,j) = G'*B(:,j);
        V(:,[i,j]) = V(:,[i,j])*G;
        if i > 1
            B(:,i-1) = R*B(:,i-1);
            U(:,i-1) = U(:,i-1)*G;
        end
    end
    if i < n
        [G,R] = planerot(B(i,i:n)',B(i+1:n,i:n));
        B(i,i:n) = B(i,i:n)*G;
        B(i+1:n,i:n) = B(i+1:n,i:n)*G;
        U(:,i+1:n) = U(:,i+1:n)*G;
        if i < m-1
            B(i:i+1,i:n) = R*B(i:i+1,i:n);
            V(:,i+1:n) = V(:,i+1:n)*G;
        end
    end
end

% SVD分解
S = diag(diag(B));
S(S == 0) = eps;
U = U(:,1:n);
V = V(:,1:n);
end

function [G,R] = planerot(x,y)
% 旋转角度计算
% 输入参数：
% x,y - 待旋转的向量
% 输出参数：
% G - 旋转矩阵
% R - 旋转后的向量

if y == 0
    G = eye(length(x));
    R = x;
elseif abs(y) > abs(x)
    t = x/y;
    s = 1/sqrt(1+t^2);
    c = s*t;
    G = [c,-s;s,c];
    R = G'*[x;y];
else
    t = y/x;
    c = 1/sqrt(1+t^2);
    s = c*t;
    G = [c,-s;s,c];
    R = G'*[x;y];
end
end

其中，`planerot`函数用于计算旋转角度，`bidiag_svd`函数用于执行双边Jacobi SVD算法。函数的输入参数为待分解的矩阵`A`，输出参数为SVD分解的结果`U`、`S`和`V`。

需要注意的是，在执行SVD分解时，对角线元素为0的情况需要进行特殊处理，将其替换为一个极小的正数`eps`。