大家好,又见面了,我是你们的朋友全栈君。如果您正在找激活码,请点击查看最新教程,关注关注公众号 “全栈程序员社区” 获取激活教程,可能之前旧版本教程已经失效.最新Idea2022.1教程亲测有效,一键激活。
Jetbrains全家桶1年46,售后保障稳定
接着LU分解继续往下,就会发展出很多相关但是并不完全一样的矩阵分解,最后对于对称正定矩阵,我们则可以给出非常有用的cholesky分解。这些分解的来源就在于矩阵本身存在的特殊的
结构。对于矩阵A,如果没有任何的特殊结构,那么可以给出A=L*U分解,其中L是下三角矩阵且对角线全部为1,U是上三角矩阵但是对角线的值任意,将U正规化成对角线为1的矩阵,产生分解A = L*D*U, D为对角矩阵。如果A为对称矩阵,那么会产生A=L*D*L分解。如果A为正定对称矩阵,那么就会产生A=G*G,可以这么理解G=L*sqrt(D)。
A=L*D*U分解对应的Matlab代码如下:
function[L, D, U] =zldu(A)
%LDU decomposition of square matrix A. The first step for Cholesky
%decomposition
[m, n] = size(A);
if m ~= n
error(‘support square matrix only’)
end
L = eye(n);
U = eye(n);
d = zeros(n,1);
for k=1:n
v = zeros(n, 1);
if k == 1
v(k:end) = A(k:end, k);
else
m = L(1:k-1, 1:k-1) \ A(1:k-1, k);
for j = 1:k-1
U(j, k) = m(j) / d(j);
end
v(k:end) = A(k:end, k) – L(k:end, 1:k-1)*m(:);
end
d(k) = v(k);
if k < n
L(k+1:end, k) = v(k+1:end)/v(k);
end
end
D = diag(d);
分解的稳定性和精度结果如下:
mean of my lu : 9.0307e-15
variance of my lu : 4.17441e-27
mean of matlab lu : 3.70519e-16
variance of matlab lu : 2.07393e-32
这里的计算是基于Gaxpy,所以稳定性和精确度相当之好。
A=L*D*L分解对应代码如下,这里要求A必须为对称矩阵:
function[D, L] =zldl(A)
%A = L*D*L’ another version of LU decomposition for matrix A
[m, n] = size(A);
if m ~= n
error(‘support square matrix only’)
end
L = eye(n);
d = zeros(n,1);
for k=1:n
v = zeros(n,1);
for j=1:k-1
v(j) = L(k, j)*d(j);
end
v(k) = A(k,k) – L(k, 1:k-1)*v(1:k-1);
d(k) = v(k);
L(k+1:end, k) = (A(k+1:end,k) – A(k+1:end, 1:k-1)*v(1:k-1)) / v(k);
end
D = diag(d);
对应分解的精确度和稳定度如下:
mean of my lu : 35.264
variance of my lu : 29011.2
mean of matlab lu : 5.88824e-16
variance of matlab lu : 8.40037e-32
使用如下的代码做测试:
n = 1500;
my_error = zeros(1, n);
sys_error = zeros(1, n);
for i = 1:n
test = gensys(5);
[zd, zl] = zldl(test);
[l, d] = ldl(test);
my_error(i) = norm(zl*zd*(zl’) – test, ‘fro’);
sys_error(i) = norm(l*d*(l’) – test, ‘fro’);
end
fprintf(‘mean of my lu : %g\n’, mean(my_error));
fprintf(‘variance of my lu : %g\n’, var(my_error));
fprintf(‘mean of matlab lu : %g\n’, mean(sys_error));
fprintf(‘variance of matlab lu : %g\n’, var(sys_error));
对于运算的精度如此之低的原因并不清楚
A=G*G’; cholesky分解对应的代码如下:
function[G] =zgaxpychol(A)
%cholesky decomposition for symmetric positive definite matrix
%the only requirement is matrix A: symmetric positive definite
[m, n] = size(A);
if m ~= n
error(‘support square matrix only’)
end
G = eye(n);
for k=1:n
v = A(:,k);
if k > 1
v(:) = v(:) – G(:,1:k-1)*G(k,1:k-1)’;
end
G(k:end, k) = v(k:end) / sqrt(v(k));
end
对应的测试结果如下
mean of my lu : 1.10711e-15
variance of my lu : 3.04741e-31
mean of matlab lu : 5.5205e-16
variance of matlab lu : 9.64928e-32
自己代码的精确度和稳定性可以媲美Matlab的代码,产生这种结果的原因应该是positive sysmetric definite matrix的原因,这段代码基于gaxpy的结果,下面给出另外一种基于外积的运算结果。
function[G] =zopchol(A)
%cholesky decomposition based on rank-1 matrix update
[m, n] = size(A);
if m ~= n
error(‘support square matrix only’)
end
G = zeros(n);
for k=1:n
G(k,k) = sqrt(A(k,k));
G(k+1:end, k) = A(k+1:end, k) / G(k,k);
%update matrix A
for j = (k+1):n
A(k+1:end,j) = A(k+1:end,j) – G(j,k)*G(k+1:end,k);
end
end
对应的测试结果如下:
mean of my lu : 9.33114e-16
variance of my lu : 1.71179e-31
mean of matlab lu : 9.92241e-16
variance of matlab lu : 1.60667e-31
对应的测试程序如下,这里使用系统自带的chol函数完成cholesky分解。
n = 1500;
my_error = zeros(1, n);
sys_error = zeros(1, n);
for i = 1:n
test = genpd(5);
[zg] = zopchol(test);
l = chol(test, ‘lower’);
my_error(i) = norm(zg*(zg’) – test, ‘fro’);
sys_error(i) = norm(l*(l’) – test, ‘fro’);
end
fprintf(‘mean of my lu : %g\n’, mean(my_error));
fprintf(‘variance of my lu : %g\n’, var(my_error));
fprintf(‘mean of matlab lu : %g\n’, mean(sys_error));
fprintf(‘variance of matlab lu : %g\n’, var(sys_error));
将两个结果想比较,可以发现两个版本的cholesky分解的精确度和稳定度差不多。
Cholesky分解的核心在于矩阵对称正定的结构,基于LU分解的再次扩展。
发布者:全栈程序员-用户IM,转载请注明出处:https://javaforall.cn/206724.html原文链接:https://javaforall.cn
【正版授权,激活自己账号】: Jetbrains全家桶Ide使用,1年售后保障,每天仅需1毛
【官方授权 正版激活】: 官方授权 正版激活 支持Jetbrains家族下所有IDE 使用个人JB账号...