STK++ 0.9.13
STK_Householder.h
Go to the documentation of this file.
1/*--------------------------------------------------------------------*/
2/* Copyright (C) 2004-2016 Serge Iovleff, Université Lille 1, Inria
3
4 This program is free software; you can redistribute it and/or modify
5 it under the terms of the GNU Lesser General Public License as
6 published by the Free Software Foundation; either version 2 of the
7 License, or (at your option) any later version.
8
9 This program is distributed in the hope that it will be useful,
10 but WITHOUT ANY WARRANTY; without even the implied warranty of
11 MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
12 GNU Lesser General Public License for more details.
13
14 You should have received a copy of the GNU Lesser General Public
15 License along with this program; if not, write to the
16 Free Software Foundation, Inc.,
17 59 Temple Place,
18 Suite 330,
19 Boston, MA 02111-1307
20 USA
21
22 Contact : S..._Dot_I..._At_stkpp_Dot_org (see copyright for ...)
23*/
24
25/*
26 * Project: Algebra
27 * Purpose: Define Householder methods for Algebra classes.
28 * Author: Serge Iovleff, S..._Dot_I..._At_stkpp_Dot_org (see copyright for ...)
29 *
30 **/
31
36#ifndef STK_HOUSEHOLDER_H
37#define STK_HOUSEHOLDER_H
38
40
41namespace STK
42{
55template <class Vector>
57{
58 // compute L^{\infty} norm of X, declare result and norm2 of X
59 Real scale = x.normInf(), v1, norm2 = 0.0;
60 if (scale) // if not 0.0
61 {
62 norm2 = (x.asDerived().sub(_R(x.begin()+1,x.lastIdx()))/=scale).norm2();
63 }
64 // check if the lower part is significative
65 if (norm2 < Arithmetic<Real>::epsilon())
66 { v1 = x.front(); x.front() = 0.0; }
67 else
68 {
69 Real s, aux1 = x.front()/scale;
70 // compute v1 = P_v X and beta of the Householder vector
71 v1 = (norm2 = sign(aux1, sqrt(aux1*aux1+norm2))) * scale;
72 // compute and save beta
73 x.front() = (s = aux1-norm2)/norm2;
74 // comp v and save it
75 x.asDerived().sub(_R(x.begin()+1,x.lastIdx())) /= s;
76 }
77 return v1;
78}
79
89template< class Lhs, class Rhs>
91{
92 // compute the product
93 Real sum = x[ v.begin()] /* *1.0 */;
94 for (int i= v.begin()+1; i<v.end(); i++) sum += x[i] * v[i];
95 // return <x,v>
96 return(sum);
97}
98
108template < class Lhs, class Rhs>
110{
111 typedef typename hidden::Traits<Lhs>::Col ColVector;
112 // get beta
113 Real beta = v.front();
114 if (beta)
115 {
116 // Multiplication of the cols by P=I+beta vv'
117 for (int j=M.beginCols(); j<M.endCols(); j++)
118 {
119 // a ref on the jth column of M
120 ColVector Mj(M.asDerived(), v.range(), j);
121 // Computation of aux=beta* <v,M^j>
122 Real aux = dotHouse( Mj, v) * beta;
123 // updating columns
124 Mj.front() += aux;
125 for(int i = v.begin()+1; i < v.end(); ++i) Mj[i] += v[i] * aux;
126 }
127 }
128}
129
141template < class Lhs, class Rhs>
143{
144 typedef typename hidden::Traits<Rhs>::Col ColVector;
145 // compute the number of iterations
146 int first = H.beginCols(), last = std::min( H.lastIdxCols(), H.lastIdxRows());
147 // get range of the first Householder vector
148 Range range_ve(last, H.lastIdxRows(), 0);
149 // iterations
150 for (int j=last; j>=first; j--)
151 {
152 // apply left Householder vector to M
153 ColVector v(H.asDerived(), range_ve, j);
155 // decrease range of the Householder vector
156 range_ve.decFirst(1);
157 }
158}
159
169template < class Lhs, class Rhs>
171{
172 // get beta
173 Real beta = v.front();
174 if (beta)
175 {
176 // Multiplication of the cols by P=I+beta vv'
177 for (int i=M.beginRows(); i<M.endRows(); i++)
178 {
179 // a ref on the ith row of M
180 typename hidden::Traits<Lhs>::Row Mi(M.asDerived(), v.range(), i);
181 // Computation of aux=beta* <v,M_i>
182 Real aux = dotHouse( Mi, v) * beta;
183 // updating column
184 Mi.front() += aux;
185 // Computation of M_i + beta <v,M_i> v = M_i + aux v'
186 for (int i=v.begin()+1; i<v.end(); i++) Mi[i] += v[i] * aux;
187 }
188 }
189}
190
201template < class TContainer2D, class Rhs>
203{
204 // compute the number of iterations
205 int first = H.beginCols(), last = std::min( H.lastIdxCols(), H.lastIdxRows());
206 // get range of the first Householder vector
207 Range range_ve(last, H.lastIdxRows());
208 // iterations
209 for (int j=last; j>=first; j--)
210 {
211 // apply left Householder vector to M
212 typename Rhs::Col v(H.asDerived(), range_ve, j);
214 // decrease range of the Householder vector
215 range_ve.decFirst(1);
216 }
217}
218
219} // namespace STK
220
221#endif /*STK_HOUSEHOLDER_H*/
In this file we define the base class for Arrays.
#define _R(first, last)
Utility macro that can be used in a similar way that first:last.
Definition STK_Range.h:53
The MultidimRegression class allows to regress a multidimensional output variable among a multivariat...
Index sub-vector region: Specialization when the size is unknown.
Definition STK_Range.h:265
Real dotHouse(ExprBase< Lhs > const &x, ExprBase< Rhs > const &v)
dot product with a Householder vector.
void applyLeftHouseholderVector(ArrayBase< Lhs > const &M, ExprBase< Rhs > const &v)
left multiplication by a Householder vector.
void applyRightHouseholderVector(ArrayBase< Lhs > const &M, ExprBase< Rhs > const &v)
right multiplication by a Householder vector.
void applyLeftHouseholderArray(ArrayBase< Lhs > const &M, ArrayBase< Rhs > const &H)
left multiplication by a Householder array.
void applyRightHouseholderArray(ArrayBase< TContainer2D > const &M, ArrayBase< Rhs > const &H)
left multiplication by a Householder ArrayXX.
Real house(ArrayBase< Vector > &x)
Compute the Householder vector v of a vector x.
Arrays::SumOp< Lhs, Rhs >::result_type sum(Lhs const &lhs, Rhs const &rhs)
convenience function for summing two arrays
Type sign(Type const &x, Type const &y=Type(1))
template sign value sign(x) * y: Type should be an integral type
Definition STK_Misc.h:53
double Real
STK fundamental type of Real values.
The namespace STK is the main domain space of the Statistical ToolKit project.
Arithmetic properties of STK fundamental types.