双边Jacobi方法求解任意矩阵的SVD的Matlab代码

以下是使用双边Jacobi方法求解任意矩阵的SVD的Matlab代码：

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

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

% Bidiagonalization using Householder reflection
for k = 1:n
    x = A(k:m,k);
    [v,beta] = house(x);
    A(k:m,k:n) = A(k:m,k:n) - beta*v*(v'*A(k:m,k:n));
    U(k:m,:) = U(k:m,:) - beta*v*(v'*U(k:m,:));
    if k < n
        x = A(k,k+1:n)';
        [v,beta] = house(x);
        A(k:m,k+1:n) = A(k:m,k+1:n) - beta*(A(k:m,k+1:n)*v)*v';
        V(k:n,:) = V(k:n,:) - beta*(V(k:n,:)*v)*v';
    end
end

% Singular value decomposition using Jacobi method
tol = 1e-12;
maxiter = 1000;
for iter = 1:maxiter
    % Diagonal pivoting
    for k = n:-1:2
        [c,s] = givens(A(k-1,k-1),A(k,k-1));
        A([k-1,k],:) = [c -s; s c]*A([k-1,k],:);
        V(:,[k-1,k]) = V(:,[k-1,k])*[c -s; s c]';
    end
    for k = 1:n-1
        [c,s] = givens(A(k,k),A(k,k+1));
        A(:,[k,k+1]) = A(:,[k,k+1])*[c -s; s c];
        U(:,[k,k+1]) = U(:,[k,k+1])*[c -s; s c]';
    end
    % Check convergence
    if norm(A(n,:)) < tol
        break;
    end
end

S = diag(A);
end

function [v,beta] = house(x)
% Householder reflection
% x: input vector
% v: Householder vector
% beta: scaling factor

n = length(x);
sigma = x(2:n)'*x(2:n);
v = [1; x(2:n)];
if sigma == 0 && x(1) >= 0
    beta = 0;
else
    if x(1) >= 0
        v(1) = -v(1);
    end
    beta = 2/(sigma + v(1)^2);
    v = v/sqrt(sigma + v(1)^2);
end
end

function [c,s] = givens(a,b)
% Givens rotation
% a,b: input scalars
% c,s: Givens parameters

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

该代码使用了Householder反射实现矩阵的双边三角化，然后使用Givens旋转实现奇异值分解。可以通过调用`bidiagonal_svd`函数来求解任意矩阵的SVD。