ProteoWizard
qr.hpp
Go to the documentation of this file.
1//
2// $Id$
3//
4// Original author: Robert Burke <robert.burke@cshs.org>
5//
6// Copyright 2006 Louis Warschaw Prostate Cancer Center
7// Cedars Sinai Medical Center, Los Angeles, California 90048
8//
9// Licensed under the Apache License, Version 2.0 (the "License");
10// you may not use this file except in compliance with the License.
11// You may obtain a copy of the License at
12//
13// http://www.apache.org/licenses/LICENSE-2.0
14//
15// Unless required by applicable law or agreed to in writing, software
16// distributed under the License is distributed on an "AS IS" BASIS,
17// WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
18// See the License for the specific language governing permissions and
19// limitations under the License.
20//
21
22#ifndef _QR_HPP
23#define _QR_HPP
24
25#include <boost/numeric/ublas/matrix.hpp>
26#include <boost/numeric/ublas/matrix_proxy.hpp>
27#include <boost/numeric/ublas/vector.hpp>
28#include <boost/numeric/ublas/io.hpp>
29
30namespace pwiz {
31namespace math {
32
33
34// Constructs a matrix to reflect a vector x onto ||x|| * e1.
35//
36// \param x vector to reflect
37// \param F matrix object to construct reflector with
38template<class matrix_type, class vector_type>
39void Reflector(const vector_type& x, matrix_type& F)
40{
41 using namespace boost::numeric::ublas;
42
43 typedef typename matrix_type::value_type value_type;
44
45 unit_vector<value_type> e1(x.size(), 0);
46
47 //v_k = -sgn( x(1) ) * inner_prod(x) * e1 + x;
48 double x_2 = norm_2(x);
49 boost::numeric::ublas::vector<value_type>
50 v_k((x(0) >= 0 ? x_2 : -1 * x_2) * e1 + x);
51
52 //v_k = v_k / norm_2(v_k);
53 double norm_vk = norm_2(v_k);
54 if (norm_vk != 0)
55 v_k /= norm_2(v_k);
56
57 // F = A(k:m,k:n) - 2 * outer_prod(v_k, v_k) * A(k:m,k:n)
58 identity_matrix<value_type> eye(v_k.size());
59 F = matrix_type(v_k.size(), v_k.size());
60
61 F = eye - 2. * outer_prod(v_k, v_k);
62}
63
64// Returns a matrix to reflect x onto ||x|| * e1.
65//
66// \param x vector to reflect
67// \return Householder reflector for x
68template<class matrix_type, class vector_type>
69matrix_type Reflector(const vector_type& x)
70{
71 using namespace boost::numeric::ublas;
72
73 matrix_type F(x.size(), x.size());
74
75 Reflector<matrix_type, vector_type>(x, F);
76
77 return F;
78}
79
80template<class matrix_type>
81void qr(const matrix_type& A, matrix_type& Q, matrix_type& R)
82{
83 using namespace boost::numeric::ublas;
84
85 typedef typename matrix_type::size_type size_type;
86 typedef typename matrix_type::value_type value_type;
87
88 // TODO resize Q and R to match the needed size.
89 int m=A.size1();
90 int n=A.size2();
91
92 identity_matrix<value_type> ident(m);
93 if (Q.size1() != ident.size1() || Q.size2() != ident.size2())
94 Q = matrix_type(m, m);
95 Q.assign(ident);
96
97 R.clear();
98 R = A;
99
100 for (size_type k=0; k< R.size1() && k<R.size2(); k++)
101 {
102 slice s1(k, 1, m - k);
103 slice s2(k, 0, m - k);
104 unit_vector<value_type> e1(m - k, 0);
105
106 // x = A(k:m, k);
107 matrix_vector_slice<matrix_type> x(R, s1, s2);
108 matrix_type F(x.size(), x.size());
109
110 Reflector(x, F);
111
112 matrix_type temp = subrange(R, k, m, k, n);
113 //F = prod(F, temp);
114 subrange(R, k, m, k, n) = prod(F, temp);
115
116 // <<---------------------------------------------->>
117 // forming Q
118 identity_matrix<value_type> iqk(A.size1());
119 matrix_type Q_k(iqk);
120
121 subrange(Q_k, Q_k.size1() - F.size1(), Q_k.size1(),
122 Q_k.size2() - F.size2(), Q_k.size2()) = F;
123
124 Q = prod(Q, Q_k);
125 }
126}
127
128}
129}
130
131#endif // _QR_HPP
F
Definition Chemistry.hpp:81
#define A
KernelTraitsBase< Kernel >::space_type::abscissa_type x
Kernel k
void Reflector(const vector_type &x, matrix_type &F)
Definition qr.hpp:39
void qr(const matrix_type &A, matrix_type &Q, matrix_type &R)
Definition qr.hpp:81