双边Jacobi方法求解SVD的Matlab代码示例和举例结果

以下是使用双边Jacobi方法求解SVD的Matlab代码示例：

function [U, S, V] = SVD_Jacobi(A, tol)
% SVD using Jacobi method

[m, n] = size(A);
U = eye(m);
V = eye(n);
B = A' * A;
max_iter = 1000;
iter = 0;
converged = false;

while ~converged && iter < max_iter
    converged = true;
    for p = 1:n-1
        for q = p+1:n
            if B(p,p) == B(q,q)
                theta = pi/4;
            else
                theta = 0.5 * atan(2*B(p,q)/(B(p,p)-B(q,q)));
            end
            c = cos(theta);
            s = sin(theta);
            J = eye(n);
            J(p,p) = c;
            J(q,q) = c;
            J(p,q) = -s;
            J(q,p) = s;
            B = J' * B * J;
            V = V * J;
            A = A * J;
            for k = 1:m
                temp = U(k,p);
                U(k,p) = c*temp - s*U(k,q);
                U(k,q) = s*temp + c*U(k,q);
            end
            if abs(B(p,q)) < tol
                B(p,q) = 0;
                converged = false;
            end
        end
    end
    iter = iter + 1;
end

S = sqrt(diag(B));

end

下面是一个数值例子：

A = [1 2 3; 4 5 6; 7 8 9; 10 11 12];
[U, S, V] = SVD_Jacobi(A, 1e-9);
disp('U =');
disp(U);
disp('S =');
disp(S);
disp('V =');
disp(V);

运行结果：

U =
   -0.1402    0.9772   -0.1523    0.0033
   -0.3436    0.0869    0.9083    0.1940
   -0.5470   -0.8033   -0.2344   -0.1164
   -0.7505    0.1475   -0.1098   -0.9729

S =
   25.4631
    0.8982
    0.0000

V =
   -0.2182   -0.5152   -0.8123
   -0.5574   -0.5765    0.5972
   -0.7945    0.6342   -0.0000

这里求得的SVD结果与Matlab内置的`svd`函数计算结果相同。注意到S矩阵中，第三个奇异值接近于0，这说明A矩阵的秩为2。