基于双边Jacobi的对任意矩阵进行svd的Matlab代码

以下是基于双边Jacobi的对任意矩阵进行svd的Matlab代码：

function [U,S,V] = bidiagonal_jacobi(A)
% Bidiagonal Jacobi SVD algorithm
% Input:  A - m x n matrix to be factorized
% Output: U - m x m orthogonal matrix
%         S - m x n diagonal matrix with singular values
%         V - n x n orthogonal matrix
[m,n] = size(A);
U = eye(m);
V = eye(n);
p = min(m,n);
for k = 1:p
    for j = k+1:n
        % Compute Givens rotation for (k,k) and (k,j) entries
        [c,s] = givens(A(k,k),A(k,j));
        % Apply Givens rotation to A and V
        A([k,j],k:n) = [c -s; s c]' * A([k,j],k:n);
        V(:,[k,j]) = V(:,[k,j]) * [c -s; s c];
    end
    for i = k+1:m
        % Compute Givens rotation for (k,k) and (i,k) entries
        [c,s] = givens(A(k,k),A(i,k));
        % Apply Givens rotation to A and U
        A(k:n,[k,i]) = A(k:n,[k,i]) * [c -s; s c];
        U(:,[k,i]) = U(:,[k,i]) * [c -s; s c];
    end
end
S = diag(diag(A)); % Extract singular values from diagonal entries
end

function [c,s] = givens(a,b)
% Compute Givens rotation for (a,b)
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

该函数的输入为一个m×n的矩阵A，输出为其SVD分解的三个矩阵U、S和V。其中，U和V均为正交矩阵，S为对角矩阵。该函数采用双边Jacobi算法对A进行分解，具体实现请参考代码注释。