rotations.h
Go to the documentation of this file.
1 /*!
2  *
3  *
4  * \brief Some operations for creating rotation matrices
5  *
6  *
7  *
8  *
9  * \author O. Krause
10  * \date 2011
11  *
12  *
13  * \par Copyright 1995-2017 Shark Development Team
14  *
15  * <BR><HR>
16  * This file is part of Shark.
17  * <http://shark-ml.org/>
18  *
19  * Shark is free software: you can redistribute it and/or modify
20  * it under the terms of the GNU Lesser General Public License as published
21  * by the Free Software Foundation, either version 3 of the License, or
22  * (at your option) any later version.
23  *
24  * Shark is distributed in the hope that it will be useful,
25  * but WITHOUT ANY WARRANTY; without even the implied warranty of
26  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
27  * GNU Lesser General Public License for more details.
28  *
29  * You should have received a copy of the GNU Lesser General Public License
30  * along with Shark. If not, see <http://www.gnu.org/licenses/>.
31  *
32  */
33
34 #ifndef SHARK_LINALG_ROTATIONS_H
35 #define SHARK_LINALG_ROTATIONS_H
36
37 #include <shark/LinAlg/Base.h>
38 #include <shark/Core/Random.h>
39 namespace shark{ namespace blas{
40 /**
41  * \ingroup shark_globals
42  *
43  * @{
44  */
45
46 //! \brief Generates a Householder reflection from a vector to use with applyHouseholderLeft/Right
47 //!
48 //! Given a Vector x=(x0,x1,...,xn), finds a reflection with the property
49 //! (c, 0,0,...0) = (I-beta v v^t)x
50 //! and v = (x0-c,x1,x2,...,xn)
51 template<class X, class R>
52 typename X::value_type createHouseholderReflection(
53  vector_expression<X, cpu_tag> const& x,
54  vector_expression<R, cpu_tag>& reflection
55 ){
56  SIZE_CHECK(x().size() != 0);
57  SIZE_CHECK(x().size() == reflection().size());
58
59  //special case for last iteration of QR etc
60  //by convention diagonal elements are > 0
61  if(x().size() == 1){
62  reflection()(0) = 1;
63  return 2;
64  }
65
66
67  typedef typename X::value_type Value;
68
69  double norm = norm_2(x);
70  if (x()(0) >= 0.0)
71  norm *= -1.0;
72
73  noalias(reflection()) = x;
74  reflection()(0) -= norm;
75  reflection() /= (x()(0) - norm);
76  //if pivot is close to 0, this is one->numericaly stable
77  //compared to the usual formula
78  Value beta = (norm - x()(0)) / norm;
79  return beta;
80 }
81 //\brief rotates a matrix using a householder reflection
82 //
83 //calculates A*(1-beta*xx^T)
84 template<class Mat, class R, class T, class Device>
86  matrix_expression<Mat, Device> & matrix,
87  vector_expression<R, Device> const& reflection,
88  T beta
89 ){
90  SIZE_CHECK(matrix().size2() == reflection().size());
91  SIZE_CHECK(reflection().size() != 0 );
92
93  //special case for last iteration of QR etc
94  if(reflection().size() == 1){
95  matrix() *= 1-beta;
96  return;
97  }
98
99  SIZE_CHECK(matrix().size2() == reflection().size());
100  //Ax
101  blas::vector<T> temp = prod(matrix,reflection);
102
103  //A -=beta*(Ax)x^T
104  noalias(matrix()) -= beta * outer_prod(temp,reflection);
105 }
106
107
108 /// \brief rotates a matrix using a householder reflection
109 ///
110 /// calculates (1-beta*xx^T)*A
111 template<class Mat, class R, class T, class Device>
113  matrix_expression<Mat, Device> & matrix,
114  vector_expression<R, Device> const& reflection,
115  T const& beta
116 ){
117
118  SIZE_CHECK(matrix().size1() == reflection().size());
119  SIZE_CHECK(reflection().size() != 0 );
120
121  //special case for last iteration of QR etc
122  if(reflection().size() == 1){
123  matrix()*=1-beta;
124  return;
125  }
126  //x^T A
127  blas::vector<T> temp = prod(trans(matrix),reflection);
128
129  //A -=beta*x(x^T A)
130  noalias(matrix()) -= beta * outer_prod(reflection,temp);
131 }
132
133 /// \brief rotates a matrix using a householder reflection
134 ///
135 /// calculates (1-beta*xx^T)*A
136 template<class Mat, class R, class T, class Device>
138  temporary_proxy<Mat> matrix,
139  vector_expression<R, Device> const& reflection,
140  T const& beta
141 ){
142  applyHouseholderOnTheLeft(static_cast<Mat&>(matrix),reflection,beta);
143 }
144
145 /// \brief Initializes a matrix such that it forms a random rotation matrix.
146 ///
147 /// The matrix needs to be quadratic and have the proper size
148 /// (e.g. call matrix::resize before).
149 ///
150 /// One common way to do this is using Gram-Schmidt-Orthogonalisation
151 /// on a matrix which is initialized with gaussian numbers. However, this is
152 /// highly unstable for big matrices.
153 ///
154 /// This algorithm is implemented from one of the algorithms presented in
155 /// Francesco Mezzadri "How to generate random matrices from the classical compact groups"
156 /// http://arxiv.org/abs/math-ph/0609050v2
157 ///
158 /// He gives two algorithms: the first one uses QR decomposition on the random gaussian
159 /// matrix and ensures that the signs of Q are correct by multiplying every column of Q
160 /// with the sign of the diagonal of R.
161 ///
162 /// We use another algorithm implemented in the paper which works similarly, but
163 /// reversed. We apply Householder rotations H_N H_{N-1}..H_1 where
164 /// H_1 is generated from a random vector on the n-dimensional unit sphere.
165 /// this requires less operations and is thus preferable. Also only half the
166 /// random numbers need to be generated
167 template< class MatrixT>
168 void randomRotationMatrix(random::rng_type& rng, matrix_container<MatrixT, cpu_tag>& matrixC){
169  MatrixT& matrix = matrixC();
170  SIZE_CHECK(matrix.size1() == matrix.size2());
171  SIZE_CHECK(matrix.size1() > 0);
172  size_t size = matrix.size1();
173  diag(matrix) = repeat(1.0,size);
174
175  RealVector v(size);
176  //we skip the first dimension as the rotation of a 1d vector is just the identity
177  for(std::size_t i = 2; i != size+1;++i){
178  //create the random vector on the unit-sphere for the i-dimensional subproblem
179  for(std::size_t j=0;j != i; ++j){
180  v(j) = random::gauss(rng);
181  }
182  subrange(v,0,i) /=norm_2(subrange(v,0,i));//project on sphere
183
184  double v0 = v(0);
185  v(0) += v0/std::abs(v0);
186
187  //compute new norm of v
188  //~ double normvSqr = 1.0+(1+v0)*(1+v0)-v0*v0;
189  double normvSqr = norm_sqr(subrange(v,0,i));
190
191  //apply the householder rotation to the i-th submatrix
192  applyHouseholderOnTheLeft(subrange(matrix,size-i,size,size-i,size),subrange(v,0,i),2.0/normvSqr);
193
194  }
195 }
196
197 //! Creates a random rotation matrix with a certain size using the random number generator rng.
198 RealMatrix randomRotationMatrix(random::rng_type& rng, size_t size){
199  RealMatrix mat(size,size);
200  randomRotationMatrix(rng, mat);
201  return mat;
202 }
203
204
205 /** @}*/
206 }}
207 #endif