我正在尝试实施一个户主三对角化例程来对小 (n < 10) hermitian matrices in matlab. This is what I have so far:
H = 1.0e-10 * [ 0.1386 + 0.0000i 0.0974 - 0.0260i 0.0434 + 0.0094i 0.0722 + 0.0670i 0.1128 + 0.1269i; 0.0974 + 0.0260i 0.0751 + 0.0000i 0.0288 + 0.0149i 0.0388 + 0.0616i 0.0557 + 0.1112i; 0.0434 - 0.0094i 0.0288 - 0.0149i 0.0146 + 0.0000i 0.0274 + 0.0164i 0.0444 + 0.0323i; 0.0722 - 0.0670i 0.0388 - 0.0616i 0.0274 - 0.0164i 0.0719 + 0.0000i 0.1216 + 0.0116i; 0.1128 - 0.1269i 0.0557 - 0.1112i 0.0444 - 0.0323i 0.1216 - 0.0116i 0.2105 + 0.0000i];
function [T, U] = tridi(H)
% Calculates a unitarily equivalent matrix T and unitary matrix U for a
% given Hermitian matrix H, such that U'*H*U = T and T is of tridiagonal
% form.
% Setting size and starting point of iteration
n = size(H,1); T = H; U = eye(n);
% Loop over the first n-2 columns of H
for k = 1 : n-1
% Householdertransformation on column-vector H(k+1:n,k)
x = T(k+1:n,k);
% Calculating secondary diagonal entry
normValue = norm(x);
% Initializing e1
e1 = zeros(length(x),1);
e1(1) = 1;
% phase of x1
phase = sign(x(1));
% Calculates normalized Householder-Vektor u
u = x + norm(x)*phase*e1;
u = u / norm(u);
% Updating T and U:
T(k+1:n,k+1:n) = (eye(n-k) - 2*(u*u'))*T(k+1:n,k+1:n)*(eye(n-k) - 2*(u*u'));
U(2:n, k+1:n) = -phase*(U(2:n, k+1:n) - 2 * (U(2:n, k+1:n) * u) * u');
% Setting secondary diagonal entry of T
T(k+1,k) = normValue;
T(k,k+1) = normValue;
% Setting appropriate row and column of T to zero
T(k+2:n,k) = 0;
T(k,k+2:n) = 0;
end
% Ensure that T is real
T = real(T);
end
但是，对于给定的数据集，当调用大约 150.000 次时，
T(k+1:n,k+1:n) = (eye(n-k) - 2*(u*u'))*T(k+1:n,k+1:n)*(eye(n-k) - 2*(u*u'))我已经尝试用矩阵向量乘法代替矩阵-矩阵乘法：
% Updating T and U:
% Calculating P*T implicitly
T(k+1:n, k+1:n) = T(k+1:n, k+1:n) - 2 * u * (u' * T(k+1:n, k+1:n));
% Calculating T*P and U*P implicitly
T(k+1:n, k+1:n) = T(k+1:n, k+1:n) - 2 * (T(k+1:n, k+1:n) * u) * u';
U(2:n, k+1:n) = -phase*(U(2:n, k+1:n) - 2 * (U(2:n, k+1:n) * u) * u');
尚未证明速度更快 - 可能是因为 matlab 矩阵乘法经过了如此充分的优化。
所以我现在正在努力降低运行时间（特别是与 Matlab 函数
hess.m我的一个想法 - 这篇文章最初的主题 - 是使用我知道的每次迭代 T 都是埃尔米特式的，这意味着只需要计算下三角部分和对角线；然后可以将下三角部分复制到上三角部分。
也可能只是将
Talfbet假设你有两个矩阵
ABAB显然我不想使用
A*B我熟悉函数
triltriu我只处理大小为 10 或更小的矩阵，因此使用稀疏矩阵可能不是正确的方法。
我想过只写一个循环，这样对角线下的所有元素都单独计算，存储在矩阵的下部，然后复制到上部，但考虑到Matlab文档坚持向量化代码而不使用循环，我真的不知道如何从不同的方式来处理它。
这里容易实现的目标是减少这里的重复计算。您在热循环中计算
I(n-k) - 2 * u * u'indices = k+1:n;
% Updating T and U:
I = eye(n - k);
two_uu = 2 * (u * u');
I_minus_two_uu = I - two_uu;
T(indices, indices) = I_minus_two_uu * T(indices, indices) * I_minus_two_uu;
U(2:n, indices) = -phase * U(2:n, indices) * I_minus_two_uu;
我使用您的代码创建了一个简短的脚本，以使用一些代码来分析结果，以确保更改不会干扰功能。由于浮点运算顺序的不同，数值不会相同，但它们是
tol = sqrt(eps); % 1.4901e-08
for i = 1:2e5
[T1, U1] = tridi(H);
[T2, U2] = tridi_faster(H);
% Ensure we get the exact same results
assert(all(ismembertol(T1, T2, tol), "all"));
assert(all(ismembertol(real(U1), real(U2), tol), "all"));
assert(all(ismembertol(imag(U1), imag(U2), tol), "all"));
end
这里有两个功能：
tridi()tridi_faster()H = 1.0e-10 * [ 0.1386 + 0.0000i 0.0974 - 0.0260i 0.0434 + 0.0094i 0.0722 + 0.0670i 0.1128 + 0.1269i; 0.0974 + 0.0260i 0.0751 + 0.0000i 0.0288 + 0.0149i 0.0388 + 0.0616i 0.0557 + 0.1112i; 0.0434 - 0.0094i 0.0288 - 0.0149i 0.0146 + 0.0000i 0.0274 + 0.0164i 0.0444 + 0.0323i; 0.0722 - 0.0670i 0.0388 - 0.0616i 0.0274 - 0.0164i 0.0719 + 0.0000i 0.1216 + 0.0116i; 0.1128 - 0.1269i 0.0557 - 0.1112i 0.0444 - 0.0323i 0.1216 - 0.0116i 0.2105 + 0.0000i];
function [T, U] = tridi(H)
% Calculates a unitarily equivalent matrix T and unitary matrix U for a
% given Hermitian matrix H, such that U'*H*U = T and T is of tridiagonal
% form.
% Setting size and starting point of iteration
n = size(H,1); T = H; U = eye(n);
% Loop over the first n-2 columns of H
for k = 1 : n-1
% Householdertransformation on column-vector H(k+1:n,k)
x = T(k+1:n,k);
% Calculating secondary diagonal entry
normValue = norm(x);
% Initializing e1
e1 = zeros(length(x),1);
e1(1) = 1;
% phase of x1
phase = sign(x(1));
% Calculates normalized Householder-Vektor u
u = x + norm(x)*phase*e1;
u = u / norm(u);
% Updating T and U:
T(k+1:n,k+1:n) = (eye(n-k) - 2*(u*u'))*T(k+1:n,k+1:n)*(eye(n-k) - 2*(u*u'));
U(2:n, k+1:n) = -phase*(U(2:n, k+1:n) - 2 * (U(2:n, k+1:n) * u) * u');
% Setting secondary diagonal entry of T
T(k+1,k) = normValue;
T(k,k+1) = normValue;
% Setting appropriate row and column of T to zero
T(k+2:n,k) = 0;
T(k,k+2:n) = 0;
end
% Ensure that T is real
T = real(T);
end
function [T, U] = tridi_faster(H)
% Calculates a unitarily equivalent matrix T and unitary matrix U for a
% given Hermitian matrix H, such that U'*H*U = T and T is of tridiagonal
% form.
% Setting size and starting point of iteration
n = size(H,1); T = H; U = eye(n);
% Loop over the first n-2 columns of H
for k = 1 : n-1
indices = k+1:n;
% Householdertransformation on column-vector H(k+1:n,k)
x = T(indices,k);
% Calculating secondary diagonal entry
normValue = norm(x);
% Initializing e1
e1 = zeros(length(x),1);
e1(1) = 1;
% phase of x1
phase = sign(x(1));
% Calculates normalized Householder-Vektor u
u = x + norm(x)*phase*e1;
u = u / norm(u);
% Updating T and U:
I = eye(n - k);
two_uu = 2 * (u * u');
I_minus_two_uu = I - two_uu;
T(indices,indices) = I_minus_two_uu * T(indices,indices) * I_minus_two_uu;
U(2:n, indices) = -phase * U(2:n, indices) * I_minus_two_uu;
% Setting secondary diagonal entry of T
T(k+1,k) = normValue;
T(k,k+1) = normValue;
% Setting appropriate row and column of T to zero
T(k+2:n,k) = 0;
T(k,k+2:n) = 0;
end
% Ensure that T is real
T = real(T);
end
tol = sqrt(eps); % 1.4901e-08
for i = 1:2e5
[T1, U1] = tridi(H);
[T2, U2] = tridi_faster(H);
% Ensure we get the exact same results
assert(all(ismembertol(T1, T2, tol), "all"));
assert(all(ismembertol(real(U1), real(U2), tol), "all"));
assert(all(ismembertol(imag(U1), imag(U2), tol), "all"));
end
使用内置 MATLAB 分析器进行 200,000 次迭代，得到以下结果。请记住，确切时间存在很大差异，但相对时间相似。
| 功能 | 通话 | 时间 |
|---|---|---|
| 原创 | 200,000 | 4.872秒 |
| 改善了 | 200,000 | 3.726秒 |
总体来说，有明显的加速，但是对于 5 秒内运行的代码，尝试优化代码的时间将远远超过代码的运行时间。
如果您确实需要实现一个数量级的加速，则需要使用不同的算法。也许有一种方法可以删除函数中的 for 循环，或者使用内置的 MATLAB 函数，但这需要了解函数实际执行的数学运算。
