ProteoWizard
HouseholderQR.hpp
Go to the documentation of this file.
1//
2// $Id$
3//
4//
5// Licensed under the Apache License, Version 2.0 (the "License");
6// you may not use this file except in compliance with the License.
7// You may obtain a copy of the License at
8//
9// http://www.apache.org/licenses/LICENSE-2.0
10//
11// Unless required by applicable law or agreed to in writing, software
12// distributed under the License is distributed on an "AS IS" BASIS,
13// WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
14// See the License for the specific language governing permissions and
15// limitations under the License.
16//
17//
18// Original author: Robert Burke <robert.burke@gmail.com>
19//
20// This code taken from the following site:
21// http://www.crystalclearsoftware.com/cgi-bin/boost_wiki/wiki.pl?Effective_UBLAS/Matrix_Inversion
22//
23
24#ifndef HOUSEHOLDERQR_HPP
25#define HOUSEHOLDERQR_HPP
26
27#include <boost/numeric/ublas/matrix.hpp>
28#include <boost/numeric/ublas/matrix_proxy.hpp>
29#include <boost/numeric/ublas/vector_proxy.hpp>
30#include <boost/numeric/ublas/io.hpp>
31#include <boost/numeric/ublas/matrix_expression.hpp>
32
33namespace ublas = boost::numeric::ublas;
34
35namespace pwiz {
36namespace math {
37
38template<class T>
39void TransposeMultiply (const ublas::vector<T>& vector,
40 ublas::matrix<T>& result,
41 size_t size)
42{
43 result.resize (size,size);
44 result.clear ();
45 for(unsigned int row=0; row< vector.size(); ++row)
46 {
47 for(unsigned int col=0; col < vector.size(); ++col)
48 result(row,col) = vector(col) * vector(row);
49
50 }
51}
52
53template<class T>
54void HouseholderCornerSubstraction (ublas::matrix<T>& LeftLarge,
55 const ublas::matrix<T>& RightSmall)
56{
57 using namespace boost::numeric::ublas;
58 using namespace std;
59 if(
60 !(
61 (LeftLarge.size1() >= RightSmall.size1())
62 && (LeftLarge.size2() >= RightSmall.size2())
63 )
64 )
65 {
66 cerr << "invalid matrix dimensions" << endl;
67 return;
68 }
69
70 size_t row_offset = LeftLarge.size2() - RightSmall.size2();
71 size_t col_offset = LeftLarge.size1() - RightSmall.size1();
72
73 for(unsigned int row = 0; row < RightSmall.size2(); ++row )
74 for(unsigned int col = 0; col < RightSmall.size1(); ++col )
75 LeftLarge(col_offset+col,row_offset+row) -= RightSmall(col,row);
76}
77
78template<class T>
79void HouseholderQR (const ublas::matrix<T>& M,
80 ublas::matrix<T>& Q,
81 ublas::matrix<T>& R)
82{
83 using namespace boost::numeric::ublas;
84 using namespace std;
85
86 if(
87 !(
88 (M.size1() == M.size2())
89 )
90 )
91 {
92 cerr << "invalid matrix dimensions" << endl;
93 return;
94 }
95 size_t size = M.size1();
96
97 // init Matrices
98 matrix<T> H, HTemp;
99 HTemp = identity_matrix<T>(size);
100 Q = identity_matrix<T>(size);
101 R = M;
102
103 // find Householder reflection matrices
104 for(unsigned int col = 0; col < size-1; ++col)
105 {
106 // create X vector
107 ublas::vector<T> RRowView = column(R,col);
108 vector_range< ublas::vector<T> > X2 (RRowView, range (col, size));
109 ublas::vector<T> X = X2;
110
111 // X -> U~
112 if(X(0) >= 0)
113 X(0) += norm_2(X);
114 else
115 X(0) += -1*norm_2(X);
116
117 HTemp.resize(X.size(),X.size(),true);
118
119 TransposeMultiply(X, HTemp, X.size());
120
121 // HTemp = the 2UUt part of H
122 HTemp *= ( 2 / inner_prod(X,X) );
123
124 // H = I - 2UUt
125 H = identity_matrix<T>(size);
127
128 // add H to Q and R
129 Q = prod(Q,H);
130 R = prod(H,R);
131 }
132}
133
134}
135}
136
137#endif // HOUSEHOLDERQR_HPP
138
H
Definition Chemistry.hpp:80
void HouseholderQR(const ublas::matrix< T > &M, ublas::matrix< T > &Q, ublas::matrix< T > &R)
void HouseholderCornerSubstraction(ublas::matrix< T > &LeftLarge, const ublas::matrix< T > &RightSmall)
void TransposeMultiply(const ublas::vector< T > &vector, ublas::matrix< T > &result, size_t size)
STL namespace.