cholesky分解_java toarray方法

cholesky分解_java toarray方法接着LU分解继续往下,就会发展出很多相关但是并不完全一样的矩阵分解,最后对于对称正定矩阵,我们则可以给出非常有用的cholesky分解。这些分解的来源就在于矩阵本身存在的特殊的结构。对于矩阵A,如果没有任何的特殊结构,那么可以给出A=L*U分解,其中L是下三角矩阵且对角线全部为1,U是上三角矩阵但是对角线的值任意,将U正规化成对角线为1的矩阵,产生分解A=L*D*U,D为对角矩阵。如果A为对…

大家好,又见面了,我是你们的朋友全栈君。如果您正在找激活码,请点击查看最新教程,关注关注公众号 “全栈程序员社区” 获取激活教程,可能之前旧版本教程已经失效.最新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账号...

(0)


相关推荐

  • java.lang.AbstractMethodError异常

    java.lang.AbstractMethodError异常前言:今天APP在客户华为Android10版本手机上闪退,报如下异常,10以下版本无任何问题,异常信息也没有自己APP报错的堆栈信息,无法直接定位到是那行代码导致的问题,苦逼了半天,梳理思路,终于解决了。分析:由于APP没有做Q版本适配,难道不兼容吗?于是安装Q版模拟器进行复现问题,直接运行安装发现没有闪退,难道是只有在华为手机上才出现问题吗?为了验证这个观点,就用客户相同的安装包进行…

  • 学术应用使用node-http-proxy集成谷歌学术「建议收藏」

    学术应用使用node-http-proxy集成谷歌学术「建议收藏」搞科研的同学肯定离不开谷歌学术,谷歌学术搜索是文献搜索下载一大利器。之前实验室开发了一款学术应用,遗留了历史问题,就是没有解决文献搜索的功能,而这个任务最后落在我的身上。我采用的方案就是集成谷歌学术,但是国内的网络环境,你懂的,自然状态下根本就访问不了谷歌学术的,你得翻墙才能访问。你不能期望使用你开发的学术应用都能翻墙访问谷歌学术(虽然搞科研的人电脑翻墙软件肯定都准备好了!),所以呢我还要给谷歌学术

  • linux生成md5文件_keil5如何生成bin文件

    linux生成md5文件_keil5如何生成bin文件MDK5如何生成bin文件MDK5在生成bin文件时经常会遇到找不到文件路径导致的生成失败,采用下面这条命令可以直接在.axf目录下生成对应的bin文件:再mdk的Target->User下面这条命令:fromelf.exe–bin-o“$L@L.bin”“#L”生成结果:…

    2022年10月20日
  • 指令周期的四个阶段_单片机指令周期与机器周期

    指令周期的四个阶段_单片机指令周期与机器周期时钟周期时钟周期也称为振荡周期,定义为时钟脉冲的倒数(可以这样来理解,时钟周期就是单片机外接晶振的倒数,例如12M的晶振,它的时间周期就是1/12us),是计算机中最基本的、最小的时间单位。在一个时钟周期内,CPU仅完成一个最基本的动作。对于某种单片机,若采用了1MHZ的时钟频率,则时钟周期为1us;若采用4MHZ的时钟频率,则时钟周期为250ns。由于时钟脉冲是计算机的基本工作脉冲…

    2022年10月13日
  • vim 设置搜索高亮_vim取消搜索后高亮持续

    vim 设置搜索高亮_vim取消搜索后高亮持续高亮有两种方法1)临时高亮vimfile文件输入:sethlsearch2)永久高亮vim~/.vimrc文件中添加sethlsearch取消搜索高亮:setnohlsearch或者简写的:noh

  • erlang如何有效地监视大量的并发连接

    erlang如何有效地监视大量的并发连接

    2021年12月30日

发表回复

您的电子邮箱地址不会被公开。

关注全栈程序员社区公众号