java 矩阵分解_计算方法(三)矩阵分解1-正交分解(QR分解)
正交分解
矩阵的正交分解又称为QR分解,是将矩阵分解为一个正交矩阵Q和一个上三角矩阵的乘积的形式。
任意实数方阵A,都能被分解为 。这里的Q为正交单位阵,即 R是一个上三角矩阵。这种分解被称为QR分解。 QR分解也有若干种算法,常见的包括Gram–Schmidt、Householder和Givens算法。 QR分解是将矩阵分解为一个正交矩阵与上三角矩阵的乘积。用一张图可以形象地表示QR分解:
为啥我们需要正交分解呢?
实际运用过程中,QR分解经常被用来解线性最小二乘问题,这个问题我们后面讲述。
提到正交分解就不得不讨论(Householder transformation Householder变换)豪斯霍尔德变换和(Schmidt orthogonalization Schmidt正交化)施密特正交化
Schmidt正交化
定理1 设A是n阶实非奇异矩阵,则存在正交矩阵Q和实非奇异上三角矩阵R使A有QR分解;且除去相差一个对角元素的绝对值(模)全等于1的对角矩阵因子外,分解是唯一的.
定理2 设A是m×n实矩阵,且其n个列向量线性无关,则A有分解A=QR,其中Q是m×n实矩阵,且满足QHTQ=E,R是n阶实非奇异上三角矩阵该分解除去相差一个对角元素的绝对值(模)全等于1的对角矩阵因子外是唯一的.用Schmidt正交化分解方法对矩阵进行QR分解时,所论矩阵必须是列满秩矩阵。
算法步骤
写出矩阵的列向量;
列向量按照Schmidt正交化正交;
得出矩阵的Q′,R′;
对R′的列向量单位化得到Q,R′的每行乘R′每列的模得푹
matlab代码
function[X,Q,R] = QRSchmidt(A,b)
%方阵的QR的Gram-Schmidt正交化分解法,并用于求解AX=b方程组[m,n]=size(A);
if m~=n
%如果不是方阵,则不满足QR分解要求
disp('不满足QR分解要求');
end
Q=zeros(m,n);
X=zeros(n,1);
R=zeros(n);
for k=1:nR(k,k)=norm(A(:,k));
if R(k,k)==0
break;
end
Q(:,k)=A(:,k)/R(k,k);
for j=k+1:n
R(k,j)=Q(:,k)'*A(:,j);
A(:,j)=A(:,j)-R(k,j)*Q(:,k);
end
if nargin==2
b=Q'* b;
X(n)=b(n)/R(n,n);
for i=n-1:-1:1
X(i)=(b(i)-sum(R(i,i+1:n).*X(i+1:n)'))/R(i,i);
end
else
X=[];
end
end
Householder变换
设A为任一n阶方阵,则必存在n阶酉矩阵Q和n阶上三角阵R,使得A=QR
设w∈Cn是一个单位向量,令
则称H是一个Householder矩阵或Householder变换。则对于任意的
存在Householder矩阵H,使得Hx=-au。其中
酉矩阵(unitary matrix)
若n阶复矩阵A满足
则称A为酉矩阵,记之为
其中,Ah是A的共轭转置
酉矩阵性质
如果A是酉矩阵
也是酉矩阵;
det(A)=1;
充分条件是它的n个列向量是两两正交的单位向量。
算法步骤
将矩阵A按列分块写成A=(α1,α2,...,αn).如果α1≠0,则可得,存在n阶householder矩阵H1使得
于是有
如果α1=0,则直接进行下一步,此时相当于取
,而a1=0.
将矩阵An-1按列分块写成An-1=(αi,α2,... ,αn-1)。如果α1≠0,则可得,存在n-1阶householder矩阵H’2使得
于是有
此时,令
则H2是n阶Householder矩阵,且使
如果α1=0,则直接进行下一步
对n-2阶矩阵继续进行类似的变换,如此下去,之多在第n-1步,我们可以找到Householder矩阵H1,H2,...,Hn-1使得
令
,则Q是酉矩阵之积,从而必有酉矩阵并且A=QR
matlab代码
function[ X,Q,R ] = QRHouseholder(A,b)
%用Householder变换将方阵A分解为正交Q与上三角矩阵R的乘积,并用于求解AX=b方程组
[n,n]=size(A);
E=eye(n);
X=zeros(n,1);
R=zeros(n);
P1=E;
for k=1:n-1
%构造w,使Pk=I-2ww'
s=-sign(A(k,k))* norm(A(k:n,k));
R(k,k)=-s;
if k==1
w=[A(1,1)+s,A(2:n,k)']';
else
w=[zeros(1,k-1),A(k,k)+s,A(k+1:n,k)']';
R(1:k-1,k)=A(1:k-1,k);
end
if norm(w)~=0
w=w/norm(w);
end
P=E-2*w*w';
A=P*A;
P1=P*P1;
R(1:n,n)=A(1:n,n);
end
Q=P1';
if nargin==2
b=P1*b;
X(n)=b(n)/R(n,n);
for i=n-1:-1:1
X(i)=(b(i)-sum(R(i,i+1:n).*X(i+1:n)'))/R(i,i);
end
else
X=[];
end
matlab自带方法
%产生一个3*3大小的魔方矩阵
A=magic(3)
[Q,R]=qr(A)
使用Eigen C++ Eigen提供了几种矩阵分解的方法
分解方式
Method
矩阵满足条件
计算速度
计算精度
PartialPivLU
partialPivLu()
Invertible
++
+
FullPivLU
fullPivLu()
None
-
+++
HouseholderQR
householderQr()
None
++
+
ColPivHouseholderQR
colPivHouseholderQr()
None
+
++
FullPivHouseholderQR
fullPivHouseholderQr()
None
-
+++
LLT
llt()
Positive definite
+++
+
LDLT
ldlt()
Positive or negative semidefinite
+++
++
其中HouseholderQR、ColPivHouseholderQR、FullPivHouseholderQR是我们目前要用到的QR分解方法
C++的QR分解代码为
#include
#include
using namespace Eigen;
using namespace std;
int main() {
Matrix3d A;
A<<1,1,1,
2,-1,-1,
2,-4,5;
HouseholderQR qr;
qr.compute(A);
MatrixXd R = qr.matrixQR().triangularView();
MatrixXd Q = qr.householderQ();
std::cout << "QR2(): HouseholderQR---------------------------------------------"<< std::endl;
std::cout << "A "<< std::endl <
std::cout <
std::cout << "R"<< std::endl <
std::cout << "Q "<< std::endl <
std::cout <
return 0;
}
输出
好了大功告成,为什么我要写计算方法的文章呢,虽然现在有很多的库和包给我们调用,但是我们也不能忘了代码的本质是为了解决复杂的数学问题,从根源上去理解一种计算方法有助于我们对自身代码的优化,比如这些方法我们可以把它写到FPGA和CUDA等并行或者分布式的计算当中,加速我们的计算方法,这比直接单机去调用这些库会超乎想象的快。
java 矩阵分解_计算方法(三)矩阵分解1-正交分解(QR分解)相关推荐
- crout分解计算例题_矩阵与数值计算(2)——矩阵三角分解LU、PALU、Cholesky三角分解、QR分解...
前言 矩阵分解是设计算法的主要技巧,通过分解可以将复杂问题转化为几个简单问题求解,通常完成这一转化任务的主要技巧就是矩阵分解.例如,我们知道上三角矩阵和下三角矩阵是容易求解的,或者对角矩阵是最理想的求 ...
- 机器学习中的矩阵向量求导(三) 矩阵向量求导之微分法
在机器学习中的矩阵向量求导(二) 矩阵向量求导之定义法中,我们讨论了定义法求解矩阵向量求导的方法,但是这个方法对于比较复杂的求导式子,中间运算会很复杂,同时排列求导出的结果也很麻烦.因此我们需要其他的 ...
- c++矩阵类_数据结构-JavaScript矩阵类的设计与实现
矩阵是线性代数课学习的重点内容之一,也是线性代数常见工具之一,在应用数学.统计分析.计算机科学.计算机图像处理级物理等多学科中均有应用.矩阵主要是指数据的行列排列的形式,由行row与列col所组成,在 ...
- 线性代数矩阵行列式_非平方矩阵的行列式| 使用Python的线性代数
线性代数矩阵行列式 Prerequisites: 先决条件: Defining a Matrix 定义矩阵 Determinant of a Matrix 矩阵的行列式 Note: Determina ...
- 向量 矩阵 张量_张量,矩阵和向量有什么区别?
向量 矩阵 张量 机器学习代数 (MACHINE LEARNING ALGEBRA) Algebra is an important element of mathematics and has a ...
- 多元函数严格凹 海塞矩阵正定_海森矩阵的应用:多元函数极值的判定
海森矩阵(Hessian Matrix),又译作黑塞矩阵.海瑟矩阵. 海塞矩阵等,是一个多元函数的二阶偏导数构成的方阵,描述 了函数的局部曲率.黑塞矩阵最早于19世纪由德国数学家 Ludwig Ott ...
- 用java制作扑克牌_阿里三面被挂,幸获内推,历经5轮终于拿到口碑offer(java研发)...
每一个互联网人心中都有一个大厂梦,百度.阿里巴巴.腾讯是很多互联网人梦寐以求的地方,小编也不例外.但是,BAT等一线互联网大厂并不是想进就能够进的,它对人才的技术能力和学历都是有一定要求的,所以除了学 ...
- python矩阵相加_【python矩阵相加怎么做,这可是证明python功能的大好机会】- 环球网校...
[摘要]今天的python实践内容是为了让大家了解python矩阵相加方法,对代码编程有个感性的认知.也好让大家能够理性选择,不要盲目跟从,选择适合自己当前阶段的学习内容,循序渐进,以兴趣自我探索为向 ...
- python矩阵转置_对python矩阵转置transpose的实例讲解
对python矩阵转置transpose的实例讲解 在读图片时,会用到这么的一段代码: image_vector_len = np.prod(image_size)#总元素大小,3*55*47 img ...
最新文章
- 自定义 DataLoader
- 马斯克:今年占全球发射质量65%,星舰5月或首次轨道试飞
- tableau可视化数据分析60讲(二十)-tableau格式设置
- 计算机多媒体教室维修登记册,多媒体教学管理制度
- Ubuntu 9.10 升级到ext4
- 面试题:各大公司Java后端开发面试题总结 已看1 背1 有用 链接有必要看看
- bzoj2257瓶子与燃料——最大公约数
- rails获取json内容
- 修复windows系统快捷方式图标变成白色的问题
- 图片中添加箭头【使用PPT实现】
- vsftpd的安装和使用
- 邮箱密码忘了怎么找回?电子邮箱密码怎么改和填写?
- mysql数据库结构导出word_Windows导出MySQL数据库表结构到Word文档-DBExportDoc V1.0 For MySQL...
- 为什么oracle打不开,oracle-Ora-01081_数据库打不开_错误解决方法
- 基于单片机的智能宠物喂食器设计
- 全网最强下载神器IDM使用教程:如何利用IDM加速下载百度网盘大文件
- SpringBoot: Redis 模拟高并发商品秒杀测试
- 错误“无法找到XXX.exe的调试信息,或者调试信息不匹配。未使用调试信息生成二进制文件“的解决方案
- python控制ppt定时_python自动化怎么操作ppt?
- 看雪学院-解密入门教学(二)笔记
热门文章
- 航天信息Aisino ZK-600+II 打印机驱动
- neatdm路径_trapcode tao插件下载-AE三维物体路径动画插件(Trapcode TAO) 2.1.2 官方最新版 - 河东下载站...
- 爬虫学习笔记num5
- 阿里专家倪超:支撑海量用户的阿里中间件技术
- 将你的QQ性别改为不男不女
- 基于单片机汽车超声波防盗系统设计(毕设课设资料)
- 努比亚Z7 mini刷机教程_recovery卡刷刷机教程
- tinyCTF 2014 tt3441810
- 护照读取_如何使用Android读取护照
- 什么是CPU高速缓存