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

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

function [U, S, V] = bidiagonal_svd(A)
% Bidiagonal SVD using Jacobi rotations
% A: input matrix
% U: left singular vectors
% S: singular values
% V: right singular vectors

[m,n] = size(A);
if m < n
    A = A';
    [m,n] = size(A);
    transpose = true;
else
    transpose = false;
end

% initialize bidiagonal matrix
B = zeros(n,n);
B(1,1) = norm(A(:,1));
u = A(:,1)/B(1,1);

% perform Householder bidiagonalization
for i = 1:n-1
    v = A(:,i+1);
    B(i,i+1) = u' * v;
    v = v - B(i,i+1) * u;
    B(i+1,i) = norm(v);
    u = v/B(i+1,i);
    B(i,i) = norm(A(:,i+1) - B(i,i+1)*u);
end
B(n,n) = norm(A(:,n) - B(:,n-1)'*A(:,n)*u);

% initialize orthogonal matrices
U = eye(m);
V = eye(n);

% perform SVD using Jacobi rotations
while true
    % check for convergence
    if max(max(abs(B - diag(diag(B))))) < eps
        break
    end
    
    % perform Jacobi rotation on B
    for i = 1:n-1
        for j = i+1:n
            if abs(B(i,j)) > eps
                theta = (B(j,j) - B(i,i))/(2*B(i,j));
                t = sign(theta)/(abs(theta) + sqrt(1+theta^2));
                c = 1/sqrt(1+t^2);
                s = c*t;
                
                G = [c s; -s c];
                B([i j],i:n) = G' * B([i j],i:n);
                B(:,[i j]) = B(:,[i j]) * G;
                V(:,[i j]) = V(:,[i j]) * G;
            end
        end
    end
    
    % perform Jacobi rotation on A
    for i = 1:m-1
        for j = i+1:m
            if abs(A(i,j)) > eps
                theta = (A(j,j) - A(i,i))/(2*A(i,j));
                t = sign(theta)/(abs(theta) + sqrt(1+theta^2));
                c = 1/sqrt(1+t^2);
                s = c*t;
                
                G = [c s; -s c];
                A([i j],i:n) = G' * A([i j],i:n);
                A(1:m,[i j]) = A(1:m,[i j]) * G;
                U(:,[i j]) = U(:,[i j]) * G;
            end
        end
    end
end

% extract singular values
S = diag(B);

% sort singular values in descending order
[S, idx] = sort(S, 'descend');
U = U(:,idx);
V = V(:,idx);

% transpose results if necessary
if transpose
    temp = U;
    U = V;
    V = temp;
end

end

该算法的基本思想是将输入矩阵进行双边Jacobi旋转，将其转化为一个上三角矩阵和一个下三角矩阵的乘积形式，然后根据乘积矩阵的奇异值分解求解出左奇异向量、奇异值和右奇异向量。该算法的时间复杂度为O(n^3)，适用于小规模矩阵的奇异值分解。