vmmlib  1.7.0
 All Classes Namespaces Functions Pages
qr_decomposition.hpp
1 /*
2  * Copyright (c) 2006-2014, Visualization and Multimedia Lab,
3  * University of Zurich <http://vmml.ifi.uzh.ch>,
4  * Eyescale Software GmbH,
5  * Blue Brain Project, EPFL
6  *
7  * This file is part of VMMLib <https://github.com/VMML/vmmlib/>
8  *
9  * Redistribution and use in source and binary forms, with or without
10  * modification, are permitted provided that the following conditions are met:
11  *
12  * Redistributions of source code must retain the above copyright notice, this
13  * list of conditions and the following disclaimer. Redistributions in binary
14  * form must reproduce the above copyright notice, this list of conditions and
15  * the following disclaimer in the documentation and/or other materials provided
16  * with the distribution. Neither the name of the Visualization and Multimedia
17  * Lab, University of Zurich nor the names of its contributors may be used to
18  * endorse or promote products derived from this software without specific prior
19  * written permission.
20  * THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS"
21  * AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
22  * IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE
23  * ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT OWNER OR CONTRIBUTORS BE
24  * LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR
25  * CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF
26  * SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS
27  * INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN
28  * CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE)
29  * ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE
30  * POSSIBILITY OF SUCH DAMAGE.
31  */
32 #ifndef __VMML__QR_DECOMPOSITION__HPP__
33 #define __VMML__QR_DECOMPOSITION__HPP__
34 
35 #include <vmmlib/matrix.hpp>
36 #include <vmmlib/vector.hpp>
37 #include <vmmlib/exception.hpp>
38 
39 #include <cmath>
40 #include <vector>
41 
42 /*
43 * QR decomposition using stabilized gram-schmidt
44 * A -> matrix to be factorized
45 * Q -> orthonormal
46 * Rn -> upper triangular
47 */
48 
49 namespace vmml
50 {
51 
52 template< size_t M, size_t N, typename T >
53 void qr_decompose_gram_schmidt(
54  const matrix< M, N, T >& A_,
55  matrix< M, M, T >& Q,
56  matrix< N, N, T >& R
57  )
58 {
59  Q = 0.0;
60  R = 0.0;
61 
62  // create a copy of A_ since we will change it in the algorithm
63  matrix< M, N, T > A( A_ );
64 
65  vector< M, T > a_column, q_column;
66 
67  #if 0
68  // for each column
69  for( size_t k = 0; k < N; ++k )
70  {
71  // compute norm of A's column k
72  A.get_column( k, a_column );
73  const T a_norm = a_column.length();
74 
75  R.at( k, k ) = a_norm;
76 
77  if ( a_norm == static_cast< T >( 0.0 ) )
78  break;
79 
80  Q.set_column( k, a_column / a_norm );
81 
82  for( size_t j = k+1; j < N; ++j )
83  {
84  Q.get_column( k, q_column );
85  A.get_column( j, a_column );
86 
87  R.at( k, j ) = a_column.dot( q_column );
88 
89  for( size_t i = 0; i < M; ++i )
90  {
91  A( i, j ) = A( i, j ) - R( k, j ) * Q( i, k );
92  }
93  }
94  }
95  #else
96  vector< M, T > v;
97  for( unsigned j = 0; j < N; ++j )
98  {
99  A.get_column( j, v );
100 
101  for( unsigned i = 0; i < j; ++i )
102  {
103  Q.get_column( i, q_column );
104  A.get_column( j, a_column );
105 
106  R( i, j ) = dot( q_column, a_column );
107  v -= q_column * R(i,j);
108  }
109  R(j,j) = v.length();
110  Q.set_column( j, v / R(j,j) );
111  }
112 
113 
114  #endif
115 
116 }
117 
118 
119 } // namespace vmml
120 
121 #endif
122