Code-Eli  0.3.6
least_squares.hpp
Go to the documentation of this file.
1 /*********************************************************************************
2 * Copyright (c) 2013 David D. Marshall <ddmarsha@calpoly.edu>
3 *
4 * All rights reserved. This program and the accompanying materials
5 * are made available under the terms of the Eclipse Public License v1.0
6 * which accompanies this distribution, and is available at
7 * http://www.eclipse.org/legal/epl-v10.html
8 *
9 * Contributors:
10 * David D. Marshall - initial code and implementation
11 ********************************************************************************/
12 
13 #ifndef eli_mutil_opt_least_squares_hpp
14 #define eli_mutil_opt_least_squares_hpp
15 
16 #include "eli/code_eli.hpp"
17 
18 namespace eli
19 {
20  namespace mutil
21  {
22  namespace opt
23  {
24  template<typename data1__, typename data2__, typename data3__>
25  void least_squares_uncon(Eigen::MatrixBase<data1__> &x, const Eigen::MatrixBase<data2__> &A, const Eigen::MatrixBase<data3__> &r)
26  {
27  unsigned int flags(0);
28  if (Eigen::MatrixBase<data2__>::ColsAtCompileTime==Eigen::Dynamic)
29  flags=Eigen::ComputeThinU | Eigen::ComputeThinV;
30  else
31  flags=Eigen::ComputeFullU | Eigen::ComputeFullV;
32 
33  Eigen::JacobiSVD<typename Eigen::MatrixBase<data2__>::PlainObject > svd(A, flags);
34  x=svd.solve(r);
35  }
36 
37  template<typename data1__, typename data2__, typename data3__, typename data4__, typename data5__>
38  void least_squares_eqcon(Eigen::MatrixBase<data1__> &x, const Eigen::MatrixBase<data2__> &A, const Eigen::MatrixBase<data3__> &b, const Eigen::MatrixBase<data4__> &B, const Eigen::MatrixBase<data5__> &d)
39  {
40  typedef Eigen::Matrix<typename Eigen::MatrixBase<data4__>::Scalar, Eigen::Dynamic, Eigen::Dynamic> B_matrixD;
41  B_matrixD Q, R, temp_B, temp_R;
42  typedef Eigen::Matrix<typename Eigen::MatrixBase<data2__>::Scalar, Eigen::Dynamic, Eigen::Dynamic> A_matrixD;
43  A_matrixD A1, A2, temp_A;
44  typedef Eigen::Matrix<typename Eigen::MatrixBase<data5__>::Scalar, Eigen::Dynamic, Eigen::Dynamic> d_matrixD;
45  d_matrixD y;
46  typedef Eigen::Matrix<typename Eigen::MatrixBase<data3__>::Scalar, Eigen::Dynamic, Eigen::Dynamic> b_matrixD;
47  b_matrixD z, rhs;
48  typedef Eigen::Matrix<typename Eigen::MatrixBase<data1__>::Scalar, Eigen::Dynamic, Eigen::Dynamic> x_matrixD;
49  x_matrixD x_temp(A.cols(), b.cols());
50  typename A_matrixD::Index p(B.rows()), n(A.cols());
51  #ifdef DEBUG
52  typename A_matrixD::Index m(b.rows());
53  #endif
54 
55  // build Q and R
56  Eigen::HouseholderQR<B_matrixD> qr(B.transpose());
57  Q=qr.householderQ();
58  temp_B=qr.matrixQR();
59  temp_R=Eigen::TriangularView<B_matrixD, Eigen::Upper>(temp_B);
60  R=temp_R.topRows(p);
61  assert((R.rows()==p) && (R.cols()==p));
62  assert((Q.rows()==n) && (Q.cols()==n));
63 
64  // build A1 and A2
65  temp_A=A*Q;
66  A1=temp_A.leftCols(p);
67  A2=temp_A.rightCols(n-p);
68 #ifdef DEBUG
69  assert((A1.rows()==m) && (A1.cols()==p));
70  assert((A2.rows()==m) && (A2.cols()==n-p));
71 #endif
72  assert(A1.cols()==p);
73  assert(A2.cols()==n-p);
74 
75  // solve for y
76  y=R.transpose().lu().solve(d);
77 
78  // setup the unconstrained optimization
79  rhs=b-A1*y;
80  least_squares_uncon(z, A2, rhs);
81 
82  // build the solution
83  x_temp.topRows(p)=y;
84  x_temp.bottomRows(n-p)=z;
85  x=Q*x_temp;
86  }
87  }
88  }
89 }
90 
91 #endif
Definition: math.hpp:20
void least_squares_eqcon(Eigen::MatrixBase< data1__ > &x, const Eigen::MatrixBase< data2__ > &A, const Eigen::MatrixBase< data3__ > &b, const Eigen::MatrixBase< data4__ > &B, const Eigen::MatrixBase< data5__ > &d)
Definition: least_squares.hpp:38
void least_squares_uncon(Eigen::MatrixBase< data1__ > &x, const Eigen::MatrixBase< data2__ > &A, const Eigen::MatrixBase< data3__ > &r)
Definition: least_squares.hpp:25