参考

[1]《矩阵计算》中译本第四版，P444~451。
[2]《矩阵计算六讲》，P134~140.
[3]NEW FAST AND ACCURATE JACOBI SVD ALGORITHM: I. doi:10.1137/050639193.
[4]NEW FAST AND ACCURATE JACOBI SVD ALGORITHM: II. doi:10.1137/05063920X

实现

按照[1]的经典Jacobi法实现

function [V] = myEigNorm(A,sigma)
[m,n] = size(A);
if m~=n
    
else
    V = eye(n);
    delta = sigma*norm(A,'fro')
    omi = off(A)
    while omi>delta
        for p=1:n-1
            for q=p+1:n
                [c,s,t]=jacobi(A(p,p),A(p,q),A(q,q));
                G = eye(n);
                G(p,p) = c;
                    G(q,q) = c;%conj(c);
                    G(p,q) = s;
                    G(q,p) = -s;%conj(s);
                    G
                    A = G'*A*G
                    V = V*G
            end
        end
        omi = off(A)
    end
end


function [c,s,t]=jacobi(alpha,beta,gama)
% 求解jacobi矩阵的三个主要值

if beta~=0
    tao = (gama -alpha)/(2*beta);
    if tao>=0
        t = 1/(tao+sqrt(1+tao^2));
    else
        t = 1/(tao-sqrt(1+tao^2));
    end
    c=1/sqrt(1+t^2);
    s=t*c;
else
    c=1;
    s=0;
    t=0;
end
end

function [o] = off(A)
%计算off函数
s = 0;
[m,n]=size(A);
for i=1:m
    for j=1:n
        if i~=j
            s=s+A(i,j)^2;
        end
    end
end
o = sqrt(s);

end

测试，测试的例子来自[2]的P140例3.4.2.

直接手算的特征值为
,,

直接调用MATLAB的eig函数的结果为

-3.96678784561050e+23
-8.10000976406272e+19
1.00000000000000e+40

Jacobi法测试代码：

clc;
clear;

A= [1e40 1e29 1e19;
    1e29 1e20 1e9;
    1e19 1e9 1]

[V] = myEigNorm(A,1e-14);
S = V'*A*V;
sigmas = [S(1,1),S(2,2),S(3,3)]

结果为

1.00000000000000e+40    
9.90000000000000e+19    
0.981818181818182

确实能够达到机器理论精度。
其他速度改进我后续再学习。