cholesky分解java代码_cholesky分解
接着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分解的再次扩展。
cholesky分解java代码_cholesky分解相关推荐
- cholesky分解java代码_Cholesky 分解(转)
Cholesky 分解是把一个对称正定的矩阵表示成一个下三角矩阵L和其转置的乘积的分解. 它要求矩阵的所有特征值必须大于零,故分解的下三角的对角元也是大于零的. Cholesky分解法又称平方根法,是 ...
- c# lu分解的代码_LU分解(1)
1/6 LU 分解 LU 分解可以写成A = LU,这里的L代表下三角矩阵,U代表上三角矩阵.对应的matlab代码如下: function[L, U] =zlu(A) % ZLU - LU deco ...
- qr分解java代码_[转载]【原创】基于Givens变换QR分解Matlab代码
如果有任何问题.建议,或者更多资源.代码.视频,欢迎您访问专业的Matlab技术交流平台--Matlab技术论坛http://www.matlabsky.com Wish my dear girl,h ...
- 递归回溯--数字分解java代码
/* * 5.数字分解 * * 6 * 6 * 5+1 * 4+2 * 4+1+1 * 3+3 * 3+2+1 * 3+1+1+1 * 2+2+2 ...
- 小波分解与小波包分解代码_分解的功能参数和代码可维护性
小波分解与小波包分解代码 Code keeps changing, there's no doubt about that. We always do our best to set some roc ...
- 面对困难的代码,分解困难各个击破
怎么说呢,最近遇到一大堆程序的时候真的有点手足无措,感觉很被动阿. 迎难而上才能有所得.这次的经验是面对一大堆要写的代码该怎么办呢? 答案:先搭框架! 其实上面这四个字着实有点太忽悠人,写代码少的人很 ...
- c语言素数筛法与分解素因数,质因数分解及代码:
计算方法 短除法 求一个数分解质因数,要从最小的质数除起,一直除到结果为质数为止.分解质因数的算式的叫短除法,和除法的性质差不多,还可以用来求多个个数的公因式: 求最大公因数的一种方法,也可用来求最小 ...
- LMD改进的局部均值分解matlab代码模版
LMD改进的局部均值分解matlab代码模版 ID:9520675130845877比邻星机灵的薯片
- 颜色分类识别代码matlab——分解RGB通道通过阈值来判断
颜色分类识别代码matlab--分解RGB通道通过阈值来判断 代码下载链接 代码下载链接 代码下载链接 首先有一张包含多种颜色的图片: 然后可以编写代码,来实现分别提取出不同颜色的操作:
- 【二】Java编程之分解质因素
质因数百度百科解释 质因数(素因数或质因子)在数论里是指能整除给定正整数的质数.除了1以外,两个没有其他共同质因子的正整数称为互质.因为1没有质因子,1与任何正整数(包括1本身)都是互质.正整数的因数 ...
最新文章
- 方法 retrun 异步的值,创建一个变量直接等于一个异步方法返回的值
- C语言编写的PHP框架--yaf入门编程
- c 多线程mysql_多线程读写mysql数据库
- Windos消息驱动
- html5怎么设置勾选,word文档怎么设置输入勾选框
- java实验七输入输出流_实验六_Java的输入输出流
- [PHP] 用JSON 传输图片源码
- 蓝桥杯2015初赛-牌型种数-dfs
- 亿佰特物联网开关电源模块:压电发声器驱动器
- 流程图伪代码计算机语言,流程图与伪代码 PPT课件
- c/c++中const用法总结
- Flask-login 原理
- Windows 适配 Apple Magic TrackPad2
- Android视频录制
- 数据结构哈夫曼树(C语言版)
- 复旦大学硕士盲审 计算机学院,《复旦大学论文抽检、盲审工作的通知.doc
- 三五族异质结的自发极化、压电极化及2DEG
- 英伟达驱动安装(华为makebook13-2018 NVIDIA-GeForce-MX150)
- 用docker-compose来搭建Hadoop(一)——创建三个ubu
- 机器学习十大算法的简单介绍