头文件：
/*
* 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 QR
Matrix<complex<Type> > QR;
// diagonal of R
Vector<complex<Type> > diagR;
// constants for generating Householder vector
Vector<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 deficient
if( !isFullRank() )
return Vector<complex<Type> >();
Vector<complex<Type> > x = b;
// compute y = Q^H * b
for( 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 x
Vector<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 deficient
if( !isFullRank() )
return Matrix<complex<Type> >(0,0);
int nx = B.cols();
Matrix<complex<Type> > X = B;
// compute Y = Q^H*B
for( 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 X
Matrix<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;
}
else
cout << " 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.
扫一扫
2490
被折叠的 条评论
为什么被折叠?
到【灌水乐园】发言
