复数矩阵QR分解算法的C++实现
2019独角兽企业重金招聘Python工程师标准>>>
头文件:
/** Copyright (c) 2008-2011 Zhang Ming (M. Zhang), zmjerry@163.com** This program is free software; you can redistribute it and/or modify it* under the terms of the GNU General Public License as published by the* Free Software Foundation, either version 2 or any later version.** Redistribution and use in source and binary forms, with or without* modification, are permitted provided that the following conditions are met:** 1. Redistributions of source code must retain the above copyright notice,* this list of conditions and the following disclaimer.** 2. Redistributions in binary form must reproduce the above copyright* notice, this list of conditions and the following disclaimer in the* documentation and/or other materials provided with the distribution.** This program is distributed in the hope that it will be useful, but WITHOUT* ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or* FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License for* more details. A copy of the GNU General Public License is available at:* http://www.fsf.org/licensing/licenses*//****************************************************************************** cqrd.h** Class template of QR decomposition for complex matrix.** For an m-by-n complex matrix A, the QR decomposition is an m-by-m unitary* matrix Q and an m-by-n upper triangular matrix R so that A = Q*R.** For economy size, denotes p = min(m,n), then Q is m-by-p, and R is n-by-p,* this file provides the economic decomposition format.** The QR decompostion always exists, even if the matrix does not have full* rank, so the constructor will never fail. The Q and R factors can be* retrived via the getQ() and getR() methods. Furthermore, a solve() method* is provided to find the least squares solution of Ax=b or AX=B using the* QR factors.** Zhang Ming, 2010-12, Xi'an Jiaotong University.*****************************************************************************/#ifndef CQRD_H
#define CQRD_H#include <matrix.h>namespace splab
{template <typename Type>class CQRD{public:CQRD();~CQRD();void dec( const Matrix<complex<Type> > &A );bool isFullRank() const;Matrix<complex<Type> > getQ();Matrix<complex<Type> > getR();Vector<complex<Type> > solve( const Vector<complex<Type> > &b );Matrix<complex<Type> > solve( const Matrix<complex<Type> > &B );private:// internal storage of QRMatrix<complex<Type> > QR;// diagonal of RVector<complex<Type> > diagR;// constants for generating Householder vectorVector<Type> betaR;};// class CQRD#include <cqrd-impl.h>}
// namespace splab#endif
// CQRD_H
实现文件:
/** Copyright (c) 2008-2011 Zhang Ming (M. Zhang), zmjerry@163.com** This program is free software; you can redistribute it and/or modify it* under the terms of the GNU General Public License as published by the* Free Software Foundation, either version 2 or any later version.** Redistribution and use in source and binary forms, with or without* modification, are permitted provided that the following conditions are met:** 1. Redistributions of source code must retain the above copyright notice,* this list of conditions and the following disclaimer.** 2. Redistributions in binary form must reproduce the above copyright* notice, this list of conditions and the following disclaimer in the* documentation and/or other materials provided with the distribution.** This program is distributed in the hope that it will be useful, but WITHOUT* ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or* FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License for* more details. A copy of the GNU General Public License is available at:* http://www.fsf.org/licensing/licenses*//****************************************************************************** qrd-impl.h** Implementation for CQRD class.** Zhang Ming, 2010-12, Xi'an Jiaotong University.*****************************************************************************//*** constructor and destructor*/
template<typename Type>
CQRD<Type>::CQRD()
{
}template<typename Type>
CQRD<Type>::~CQRD()
{
}/*** Create a QR factorization for a complex matrix A.*/
template <typename Type>
void CQRD<Type>::dec( const Matrix<complex<Type> > &A )
{int m = A.rows(),n = A.cols(),p = min(m,n);Type absV0, normV, beta;complex<Type> alpha;diagR = Vector<complex<Type> >(p);betaR = Vector<Type>(p);QR = A;// main loop.for( int k=0; k<p; ++k ){// Form k-th Householder vector.normV = 0;for( int i=k; i<m; ++i )normV += norm(QR[i][k]);normV = sqrt(normV);absV0 = abs(QR[k][k]);alpha = -normV * QR[k][k]/absV0;beta = 1 / ( normV * (normV+absV0) );QR[k][k] -= alpha;// Apply transformation to remaining columns.for( int j=k+1; j<n; ++j ){complex<Type> s = 0;for( int i=k; i<m; ++i )s += conj(QR[i][k]) * QR[i][j];s *= beta;for( int i=k; i<m; ++i )QR[i][j] -= s*QR[i][k];}diagR[k] = alpha;betaR[k] = beta;}// int m = A.rows(),
// n = A.cols(),
// p = min(m,n);
//
// Type absV0, normV, beta;
// complex<Type> alpha;
//
// diagR = Vector<complex<Type> >(p);
// betaR = Vector<Type>(p);
// QR = A;
//
// // main loop.
// for( int k=0; k<p; ++k )
// {
// // Form k-th Householder vector.
// Vector<complex<Type> > v(m-k);
// for( int i=k; i<m; ++i )
// v[i-k] = QR[i][k];
//
// absV0 = abs(v[0]);
// normV = norm(v);
// alpha = -normV * v[0]/absV0;
// beta = 1 / ( normV * (normV+absV0) );
// v[0] -= alpha;
//
// for( int i=k; i<m; ++i )
// QR[i][k] = v[i-k];
//
// // Apply transformation to remaining columns.
// for( int j=k+1; j<n; ++j )
// {
// complex<Type> s = 0;
// for( int i=k; i<m; ++i )
// s += conj(v[i-k]) * QR[i][j];
// for( int i=k; i<m; ++i )
// QR[i][j] -= beta*s*v[i-k];
// }
//
// diagR[k] = alpha;
// betaR[k] = beta;
// }
}/*** Flag to denote the matrix is of full rank.*/
template <typename Type>
inline bool CQRD<Type>::isFullRank() const
{for( int j=0; j<diagR.dim(); ++j )if( abs(diagR[j]) == Type(0) )return false;return true;
}/*** Return the upper triangular factorof the QR factorization.*/
template <typename Type>
Matrix<complex<Type> > CQRD<Type>::getQ()
{int m = QR.rows(),p = betaR.dim();Matrix<complex<Type> >Q( m, p );
// for( int i=0; i<p; ++i )
// Q[i][i] = 1;for( int k=p-1; k>=0; --k ){Q[k][k] = 1 - betaR[k]*norm(QR[k][k]);for( int i=k+1; i<m; ++i )Q[i][k] = -betaR[k]*QR[i][k]*conj(QR[k][k]);for( int j=k+1; j<p; ++j ){complex<Type> s = 0;for( int i=k; i<m; ++i )s += conj(QR[i][k])*Q[i][j];s *= betaR[k];for( int i=k; i<m; ++i )Q[i][j] -= s * QR[i][k];}}return Q;
}/*** Return the orthogonal factorof the QR factorization.*/
template <typename Type>
Matrix<complex<Type> > CQRD<Type>::getR()
{int n = QR.cols(),p = diagR.dim();Matrix<complex<Type> > R( p, n );for( int i=0; i<p; ++i ){R[i][i] = diagR[i];for( int j=i+1; j<n; ++j )R[i][j] = QR[i][j];}return R;
}/*** Least squares solution of A*x = b, where A and b are complex.* Return x: a n-length vector that minimizes the two norm* of Q*R*X-B. If B is non-conformant, or if QR.isFullRank()* is false, the routinereturns a null (0-length) vector.*/
template <typename Type>
Vector<complex<Type> > CQRD<Type>::solve( const Vector<complex<Type> > &b )
{int m = QR.rows(),n = QR.cols();assert( b.dim() == m );// matrix is rank deficientif( !isFullRank() )return Vector<complex<Type> >();Vector<complex<Type> > x = b;// compute y = Q^H * bfor( int k=0; k<n; ++k ){complex<Type> s = 0;for( int i=k; i<m; ++i )s += conj(QR[i][k])*x[i];s *= betaR[k];for( int i=k; i<m; ++i )x[i] -= s * QR[i][k];}// solve R*x = y;for( int k=n-1; k>=0; --k ){x[k] /= diagR[k];for( int i=0; i<k; ++i )x[i] -= x[k]*QR[i][k];}// return n portion of xVector<complex<Type> > x_(n);for( int i=0; i<n; ++i )x_[i] = x[i];return x_;
}/*** Least squares solution of A*X = B, where A and B are complex.* return X: a matrix that minimizes the two norm of Q*R*X-B.* If B is non-conformant, or if QR.isFullRank() is false, the* routinereturns a null (0) matrix.*/
template <typename Type>
Matrix<complex<Type> > CQRD<Type>::solve( const Matrix<complex<Type> > &B )
{int m = QR.rows();int n = QR.cols();assert( B.rows() == m );// matrix is rank deficientif( !isFullRank() )return Matrix<complex<Type> >(0,0);int nx = B.cols();Matrix<complex<Type> > X = B;// compute Y = Q^H*Bfor( int k=0; k<n; ++k )for( int j=0; j<nx; ++j ){complex<Type> s = 0;for( int i=k; i<m; ++i )s += conj(QR[i][k])*X[i][j];s *= betaR[k];for( int i=k; i<m; ++i )X[i][j] -= s*QR[i][k];}// solve R*X = Y;for( int k=n-1; k>=0; --k ){for( int j=0; j<nx; ++j )X[k][j] /= diagR[k];for( int i=0; i<k; ++i )for( int j=0; j<nx; ++j )X[i][j] -= X[k][j]*QR[i][k];}// return n x nx portion of XMatrix<complex<Type> > X_( n, nx );for( int i=0; i<n; ++i )for( int j=0; j<nx; ++j )X_[i][j] = X[i][j];return X_;
}
测试代码:
/****************************************************************************** cqrd_test.cpp** CQRD class testing.** Zhang Ming, 2010-12, Xi'an Jiaotong University.*****************************************************************************/#define BOUNDS_CHECK#include <iostream>
#include <iomanip>
#include <cqrd.h>using namespace std;
using namespace splab;typedef double Type;
const int M = 4;
const int N = 3;int main()
{Matrix<Type> A(M,N);A[0][0] = 1; A[0][1] = 2; A[0][2] = 3;A[1][0] = 1; A[1][1] = 2; A[1][2] = 4;A[2][0] = 1; A[2][1] = 9; A[2][2] = 7;A[3][0] = 5; A[3][1] = 6; A[3][2] = 8;
// A = trT(A);Matrix<complex<Type> > cA = complexMatrix( A, elemMult(A,A) );CQRD<Type> qr;qr.dec(cA);Matrix<complex<Type> > Q = qr.getQ();Matrix<complex<Type> > R = qr.getR();cout << setiosflags(ios::fixed) << setprecision(3);cout << "The original matrix cA : " << cA << endl;cout << "The orthogonal matrix Q : " << Q << endl;cout << "The upper triangular matrix R : " << R << endl;cout << "Q^H * Q : " << trMult(Q,Q) << endl;cout << "cA - Q*R : " << cA - Q*R << endl;Vector<Type> b(M);b[0]= 1; b[1] = 0; b[2] = 1, b[3] = 2;Vector<complex<Type> > cb = complexVector(b);if( qr.isFullRank() ){Vector<complex<Type> > x = qr.solve(cb);cout << "The constant vector cb : " << cb << endl;cout << "The least squares solution of cA * x = cb : " << x << endl;Matrix<complex<Type> > X = qr.solve( eye( M, complex<Type>(1) ) );cout << "The least squares solution of cA * X = I : " << X << endl;cout << "The cA * X: " << cA*X << endl;}elsecout << " The matrix is rank deficient! " << endl;return 0;
}
运行结果:
The original matrix cA : size: 4 by 3
(1.000,1.000) (2.000,4.000) (3.000,9.000)
(1.000,1.000) (2.000,4.000) (4.000,16.000)
(1.000,1.000) (9.000,81.000) (7.000,49.000)
(5.000,25.000) (6.000,36.000) (8.000,64.000)The orthogonal matrix Q : size: 4 by 3
(-0.055,0.000) (-0.029,-0.000) (0.370,-0.040)
(-0.055,0.000) (-0.029,-0.000) (0.913,-0.143)
(-0.055,0.000) (-0.985,-0.158) (-0.042,-0.000)
(-0.828,-0.552) (0.043,0.039) (-0.063,-0.030)The upper triangular matrix R : size: 3 by 3
(-18.111,-18.111) (-25.565,-31.418) (-42.737,-52.676)
(0.000,0.000) (-20.071,-77.269) (-11.956,-45.432)
(0.000,0.000) (0.000,0.000) (-0.593,12.784)Q^H * Q : size: 3 by 3
(1.000,0.000) (0.000,0.000) (-0.000,-0.000)
(0.000,-0.000) (1.000,0.000) (-0.000,0.000)
(-0.000,0.000) (-0.000,-0.000) (1.000,0.000)cA - Q*R : size: 4 by 3
(-0.000,-0.000) (-0.000,-0.000) (-0.000,-0.000)
(0.000,0.000) (-0.000,-0.000) (-0.000,-0.000)
(0.000,0.000) (0.000,0.000) (0.000,-0.000)
(0.000,0.000) (0.000,0.000) (0.000,-0.000)The constant vector cb : size: 4 by 1
(1.000,0.000)
(0.000,0.000)
(1.000,0.000)
(2.000,0.000)The least squares solution of cA * x = cb : size: 3 by 1
(-0.002,-0.035)
(-0.002,-0.002)
(0.007,-0.016)The least squares solution of cA * X = I : size: 3 by 4
(-0.007,0.048) (-0.025,0.120) (-0.002,0.012) (0.004,-0.048)
(-0.001,0.017) (-0.004,0.042) (0.001,-0.014) (-0.001,-0.002)
(0.002,-0.029) (0.008,-0.072) (0.000,0.003) (0.003,0.005)The cA * X: size: 4 by 4
(0.143,-0.000) (0.348,0.017) (0.016,-0.003) (0.022,-0.016)
(0.348,-0.017) (0.858,0.000) (-0.007,0.001) (-0.009,0.007)
(0.016,0.003) (-0.007,-0.001) (1.000,-0.000) (-0.000,0.000)
(0.022,0.016) (-0.009,-0.007) (-0.000,-0.000) (0.999,0.000)Process returned 0 (0x0) execution time : 0.109 s
Press any key to continue.
转载于:https://my.oschina.net/zmjerry/blog/3757
复数矩阵QR分解算法的C++实现相关推荐
- 复数矩阵Cholesky分解算法的C++实现
2019独角兽企业重金招聘Python工程师标准>>> 头文件: /** Copyright (c) 2008-2011 Zhang Ming (M. Zhang), zmjerry ...
- Julia 矩阵QR分解和特征值
Julia 矩阵QR分解和特征值 前言 1. 施密特正交 (1) 利用施密特正交求出正交矩阵Q (2) 求出上三角矩阵R (3) 改进的消减QR分解 2. 完全QR分解 3. 矩阵QR分解的作用 (1 ...
- C语言实现实数和复数矩阵及其各种运算(四)
一.前言 本章开始,笔者将详细讲解实数和复数矩阵的模(范数).协方差等运算,并给出函数和测试demo的C代码,以及与matlab计算的结果: 并且,本章也出现了我在编写函数.后期与matlab联调跑数 ...
- 复数矩阵分解的拆解思路(矩阵求逆/特征值分解)
作者:桂. 时间:2017-10-26 07:11:02 链接:http://www.cnblogs.com/xingshansi/p/7735016.html 前言 主要记录特征值分解的硬件实现思 ...
- Python+numpy实现矩阵QR分解
感谢广东东软学院计算机系赵晨杰老师的交流. 如果实(复)非奇异矩阵A能够化成正交(酉)矩阵Q与实(复)非奇异上三角矩阵R的乘积,即A=QR,则称其为A的QR分解. Python扩展库numpy实现了矩 ...
- 将两个实数矩阵合并为一个复数矩阵
问题描述:有时需要把两个实数矩阵,一个作为实部,一个作为虚部,合并为一个复数矩阵,该如何操作? 解决办法: 假如是在第二个维度上进行合并(real: Data[:, 0, :, :] imag: Da ...
- C语言实现实数和复数矩阵及其各种运算(二)
一.前言 由于实数矩阵的运算较简单,因此在本章中,我只给出复数矩阵的相关运算,一般的实数矩阵,类似炮制即可: 复数矩阵的加/减/乘运算涉及到其复数元胞(cell)的相加减运算,由于complex.h头 ...
- C语言实现实数和复数矩阵及其各种运算(一)
一.前言 本连载文章主要内容是实现复数矩阵的各种运算,并利用matlab进行联写.联调,验证C语言编写的各个矩阵运算函数的正确性.其难点笔者认为在于矩阵的各种运算,C++中有Eigen库可用,以前在学 ...
- 复数矩阵的转置、共轭、共轭转置
复数矩阵的转置分[共轭转置]和[点转置]两种,如果不加说明默认指的是[共轭转置]. 在 MATLAB 中,点转置运算.共轭运算和共轭转置运算分别如下所示:
最新文章
- 撸了个低代码开发平台,爽!
- CCNA-EiGrp学习
- 记录运行gpu错误及解决方案
- split函数python 未定义_Python字符串方法split()中的一道坑
- Java读取模板文件您好,RtfTemplate 读取word模板生成文件
- 非常精简的Linux线程池实现(一)——使用互斥锁和条件变量
- 前端学习(1904)vue之电商管理系统电商系统之修改用户的操作
- 很用心的为你写了 9 道 MySQL 面试题,建议收藏!
- mysql 不完全插入_MySql insert插入操作不完全指北_MySQL
- 【LeetCode】【HOT】39. 组合总和(回溯)
- 为什么学习web前端,必须掌握JavaScript这门编程语言
- 修改表名的sql语句_SQL第一关——入门
- 配置ftp方式的yum源的各种排错
- Cursor finalized without prior close()
- c++按行读取txt文件中的内容,并按特定字符分割
- 微信打电话和直接打电话有什么区别吗?为什么?
- STM32 通用 Bootloader
- k3s 快速入门 - traefix 使用 - 1
- python接收163邮件以及下载附件(以163邮箱为例)
- 无人车系统(八):Udacity ‘s无人驾驶仿真环境(python与c++数据接口)
热门文章
- maven依赖avro_如何使用maven进行avro序列化
- was这么做的负载均衡_中间件(WAS、WMQ)运维 9个常见难点解析
- Windows 命令行大全
- godaddy 管理mysql_在godaddy上使用MySQL和Entity Framework的安全例...
- LeetCode - 4. 寻找两个正序数组的中位数
- 百练OJ:2767:简单密码
- Android中提示:Service Intent must be explicit:Intent
- Winform中DevExpress的TreeList的入门使用教程(附源码下载)
- 【NLP实战】Task1 数据集探索
- 强大的代码扫描工具SonarLint之安装使用