文档章节

LU分解

SVD
 SVD
发布于 2015/11/04 15:11
字数 682
阅读 246
收藏 0
function [L,U,p] = lutxloops(A)
%LU Triangular factorization
%   [L,U,p] = lup(A) produces a unit lower triangular matrix L,
%   an upper triangular matrix U and a permutation vector p,
%   so that L*U = A(p,:)

[n,n] = size(A);
p = (1:n)';

for k = 1:n-1

   % Find index of largest element below diagonal in k-th column
   m = k;
   for i = k+1:n
      if abs(A(i,k)) > abs(A(m,k))
         m = i;
      end
   end

   % Skip elimination if column is zero
   if (A(m,k) ~= 0)
   
      % Swap pivot row
      if (m ~= k)
         for j = 1:n;
            A([k m],j) = A([m k],j);
         end
         p([k m]) = p([m k]);
      end

      % Compute multipliers
      for i = k+1:n;
         A(i,k) = A(i,k)/A(k,k);
      end

      % Update the remainder of the matrix
      for j = k+1:n;
         for i = k+1:n;
            A(i,j) = A(i,j) - A(i,k)*A(k,j); 
         end
      end
   end
end
    
% Separate result
L = tril(A,-1) + eye(n,n);
U = triu(A);

之前分析过gauss elimination,其核心是LU分解。本段代码保存为lutxloops.m以上代码不再带有线性方程的右端项而已,本质上无区别。

L = tril(X,k) returns the elements on and below the kth diagonal of X. k = 0 is the main diagonal, k > 0 is above the main diagonal, and k < 0 is below the main diagonal.

U = triu(X) returns the upper triangular part of X.

以下代码是lutx.m,可以看到,少用了许多for循环。接下来就是测试这两个功能相同的程序,其各自的运行效率了。

function [L,U,p] = lutx(A)
%LUTX  Triangular factorization, textbook version
%   [L,U,p] = lutx(A) produces a unit lower triangular matrix L,
%   an upper triangular matrix U, and a permutation vector p,
%   so that L*U = A(p,:)



[n,n] = size(A);
p = (1:n)';

for k = 1:n-1

   % Find index of largest element below diagonal in k-th column
   [r,m] = max(abs(A(k:n,k)));
   m = m+k-1;

   % Skip elimination if column is zero
   if (A(m,k) ~= 0)
   
      % Swap pivot row
      if (m ~= k)
         A([k m],:) = A([m k],:);
         p([k m]) = p([m k]);
      end

      % Compute multipliers
      i = k+1:n;
      A(i,k) = A(i,k)/A(k,k);

      % Update the remainder of the matrix
      j = k+1:n;
      A(i,j) = A(i,j) - A(i,k)*A(k,j); 
   end
end

% Separate result
L = tril(A,-1) + eye(n,n);
U = triu(A);

以下是测试结果:

>> n = 525, A = randn(n,n); tic, lutxloops(A); toc


n =


   525


Elapsed time is 2.314373 seconds.

>> n = 525, A = randn(n,n); tic, lutx(A); toc


n =


   525


Elapsed time is 0.584623 seconds.

可见,运行效率差距是明显的。

>> n = 1920, A = randn(n,n); tic, lu(A); toc


n =


        1920


Elapsed time is 0.191691 seconds.

内置的lu函数的效率最高。这是因为内部已经将其编译成了可执行的机器码。由于matlab是解释执行的编程语言,且其自身的特点是对于矩阵和向量的运算能力较强。因此在Matlab 中应尽量用向量、矩阵运算,不要用 For 循环。

MATLAB (matrix laboratory) is a multi-paradigm numerical computing environment and fourth-generation programming language. A proprietary programming language developed by MathWorks, MATLAB allows matrix manipulations, plotting of functions and data, implementation of algorithms, creation of user interfaces, and interfacing with programs written in other languages, including C, C++, Java, Fortran and Python.

© 著作权归作者所有

共有 人打赏支持
SVD

SVD

粉丝 33
博文 207
码字总数 102940
作品 0
海淀
私信 提问
Armadillo之LU分解(LU factorisation/LU decomposition)

在armadillo库中,矩阵的LU分解(LU factorisation or LU decomposition)使用lu函数,lu函数有两个版本 1 lu(L,U,P,X) 其中X是欲进行分解的矩阵,分解生成L,U,P满足 1)P是一个置换矩阵(p...

桑梓狼狼
2014/08/01
0
0
使用javascript实现矩阵LU分解

版权声明:本文为博主原创文章,欢迎读者转载学习,转载请注明文章出处。 https://blog.csdn.net/qq_37338983/article/details/79918254 在线性代数中,LU分解是将一个矩阵分解为 L(单位下三...

Bruce_wjh
2018/04/12
0
0
各位大大,怎样在Eigen3.3.3库调用其支持的模块。

如何通过Eigen库,调用其支持的第三方模块。比如我想在Eigen中使用PaSTiX库的LU分解来代替Eigen中自带的LU分解,该如何进行呢? 1、PaSTix官网:http://pastix.gforge.inria.fr/files/READM...

huhuhaha
2017/05/24
96
0
你知道MATLAB ,但你知道 NMATH吗?

NMath是一个.NET的数学库,包含了NET平台上的面向对象数字计算的基础类。 产品特点如下: 单精度和双精度复数类 为以下四种数据类型提供全功能的向量和矩阵类:单精度浮点数,双精度浮点数,...

愤怒的小吉
2014/10/24
612
0
Eigen 3.2.0-beta1 发布,线性算术的C++模板库

这个beta版本引入了内置的稀疏矩阵,和真正的QZ分解和广义特征求解稠密矩阵的LU和QR因子分解,以及Ref<>参考类。同时修复了一些bug。 Eigen 是一个线性算术的C++模板库,包括:vectors, matr...

zino
2013/03/08
923
6

没有更多内容

加载失败,请刷新页面

加载更多

JeeSite4.x 消息管理、消息推送、消息提醒

实现统一的消息推送接口,包含PC消息、短信消息、邮件消息、微信消息等,无需让所有开发者了解消息是怎么发送出去的,只需了解消息发送接口即可。 所有推送消息均通过 MsgPushUtils 工具类发...

ThinkGem
13分钟前
3
0
OpenML

https://www.openml.org/search?type=data

shengjuntu
15分钟前
2
0
java强引用,软引用,弱引用和虚引用

先来简要说一下这四种引用的特性: 强引用:如果一个对象具有强引用,那垃圾回收器绝不会回收它 软引用:如果一个对象只具有软引用,则内存空间足够,垃圾回收器就不会回收它 弱引用:在垃圾...

woshixin
19分钟前
1
0
dubbo服务在kubernetes中对外暴露服务

一些场景下,在开发时可能存在dubbo consumer需要访问k8s中部署的dubbo provider,尤其是对于自建的kubernetes集群环境,tcp的端口很难代理,这就导致了开发发链接集群内的服务比较麻烦,这里...

上官胡闹
29分钟前
1
0
Java获取小程序带参二维码(太阳码)

获取小程序码 官方API地址 : https://developers.weixin.qq.com/miniprogram/dev/api/qrcode.html 首先使用官方提供的 接口B:适用于需要的码数量极多的业务场景 https://api.weixin.qq.com/...

回忆是如此忧伤
36分钟前
1
0

没有更多内容

加载失败,请刷新页面

加载更多

返回顶部
顶部