Lesson N.11 of July 31, 2014¶
Some experiments using Eigen3
File example17.cc
:open:
1/*
2 Some experiments using Eigen3
3
4 -I/usr/local/include/eigen3 to add search of headers for eigen3 library
5
6*/
7
8#include <iostream>
9#include <iomanip>
10#include <vector>
11#include <cmath>
12#include <string> // C++ string
13#include <exception> // std::exception
14
15// include headers for the Dense linear algebra of Eigen library
16#include <Eigen/Dense>
17
18using namespace std ;
19//using namespace Eigen ;
20
21int
22main() {
23
24 try {
25 // define a matrix od double (d) of size 2x2
26 // Matrix2d = Matrix<double,2,2>
27 Eigen::Matrix2d a ;
28 a << 1, 2,
29 3, 4 ;
30
31 // MatrixXd = Matrix<double,Dynamic,Dynamic>
32 Eigen::MatrixXd b(2,2) ;
33 b << 2, 3,
34 1, 4;
35
36 cout << "a + b =\n" << a + b << '\n' ;
37 cout << "a - b =\n" << a - b << '\n' ;
38 cout << "Doing a += b;" << '\n' ;
39 a += b;
40 cout << "Now a =\n" << a << '\n' ;
41 // Vector3d = Matrix<double,3,1>
42 Eigen::Vector3d v(1,2,3); // v << 1 << 2 << 3 ;
43 Eigen::Vector3d w(1,0,0); // v << 1 << 0 << 0 ;
44 cout << "-v + w - v =\n" << -v + w - v << '\n' ;
45 }
46 catch ( exception const & error ) {
47 cerr << "Exception number " << error.what() << "\n" ;
48 }
49 catch ( ... ) {
50 cerr << "Exception found: Unknown error\n" ;
51 }
52 cout << "END OF PROGRAM\n" ;
53
54 return 0 ;
55}
File example18.cc
:open:
1/*
2 Some experiments using Eigen3
3
4 -I/usr/local/include/eigen3 to add search of headers for eigen3 library
5
6
7 Linear regression, bad way and good way
8
9*/
10
11#include <iostream>
12#include <iomanip>
13#include <vector>
14#include <cmath>
15#include <string> // C++ string
16#include <exception> // std::exception
17
18// include headers for the Dense linear algebra of Eigen library
19#include <Eigen/Dense>
20
21using namespace std ;
22//using namespace Eigen ;
23
24template <typename T>
25T
26func( T x )
27{ return 10*x-3 ; }
28
29template <typename T>
30void
31generate_table( vector<T> & x,
32 vector<T> & y,
33 int N ) {
34 x.clear() ;
35 y.clear() ;
36 x.reserve(N) ;
37 y.reserve(N) ;
38
39 for ( int i = 0 ; i < N ; ++i ) {
40 T xx = ((rand() % 10000)-5000) * 0.001 ;
41 T yy = func(xx) ;
42 yy += ((rand() % 10000)-5000) * 0.0001 ;
43 x.push_back(xx) ;
44 y.push_back(yy) ;
45 }
46
47}
48
49typedef double valueType ;
50
51int
52main() {
53
54 try {
55 int N = 1000000 ;
56 vector<valueType> x, y ;
57 generate_table( x, y, N ) ;
58
59 // define a matrix od double (d) of size 2x2
60 // Matrix2d = Matrix<double,2,2>
61 Eigen::Matrix2d A ;
62 Eigen::Vector2d b, sol ;
63
64 // sum xi*xi
65 valueType bf = 0 ;
66 for ( int i = 0 ; i < N ; ++i ) bf += x[i]*x[i] ;
67 A(0,0) = bf ;
68
69 // sum xi
70 bf = 0 ;
71 for ( int i = 0 ; i < N ; ++i ) bf += x[i] ;
72 A(1,0) = A(0,1) = bf ;
73 A(1,1) = N ;
74
75 // sum xi*yi
76 bf = 0 ;
77 for ( int i = 0 ; i < N ; ++i ) bf += x[i]*y[i] ;
78 b(0) = bf ;
79
80 // sum yi
81 bf = 0 ;
82 for ( int i = 0 ; i < N ; ++i ) bf += y[i] ;
83 b(1) = bf ;
84
85 // use Cholewsky decomposition
86 Eigen::LLT<Eigen::Matrix2d> lltOfA(A) ;
87 sol = lltOfA.solve(b) ;
88 cout << "Solution = " << sol << '\n' ;
89
90 // solve the problem using least squares
91 Eigen::Matrix<double,Eigen::Dynamic,Eigen::Dynamic> Mat(N,2) ;
92 Eigen::Matrix<double,Eigen::Dynamic,1> Rhs(N) ;
93 for ( int i = 0 ; i < N ; ++i ) {
94 Mat(i,0) = x[i] ;
95 Mat(i,1) = 1 ;
96 Rhs(i) = y[i] ;
97 }
98
99 cout << "Mat^T * Mat =\n" << Mat.transpose() * Mat << '\n' ;
100 cout << "A =\n" << A << '\n' ;
101
102 sol = Mat.colPivHouseholderQr().solve(Rhs) ;
103
104 cout << "Solution (LS) = " << sol << '\n' ;
105
106 }
107 catch ( exception const & error ) {
108 cerr << "Exception number " << error.what() << "\n" ;
109 }
110 catch ( ... ) {
111 cerr << "Exception found: Unknown error\n" ;
112 }
113 cout << "END OF PROGRAM\n" ;
114
115 return 0 ;
116}
File example19.cc
:open:
1/*
2 Some experiments using Eigen3
3
4 -I/usr/local/include/eigen3 to add search of headers for eigen3 library
5
6
7 Linear regression, bad way and good way.
8 Using more Eigen3
9
10*/
11
12#include <iostream>
13#include <iomanip>
14#include <vector>
15#include <cmath>
16#include <string> // C++ string
17#include <exception> // std::exception
18
19// include headers for the Dense linear algebra of Eigen library
20#include <Eigen/Core>
21#include <Eigen/Dense>
22
23using namespace std ;
24//using namespace Eigen ;
25
26template <typename T>
27T
28func( T x )
29{ return 10*x-3 ; }
30
31template <typename T>
32void
33generate_table( Eigen::Matrix<T,Eigen::Dynamic,1> & x,
34 Eigen::Matrix<T,Eigen::Dynamic,1> & y,
35 int N ) {
36 x.resize(N) ;
37 y.resize(N) ;
38
39 for ( int i = 0 ; i < N ; ++i ) {
40 T xx = ((rand() % 10000)-5000) * 0.001 ;
41 T yy = func(xx) ;
42 yy += ((rand() % 10000)-5000) * 0.0001 ;
43 x(i) = xx ;
44 y(i) = yy ;
45 }
46
47}
48
49typedef double valueType ;
50
51int
52main() {
53
54 try {
55 int N = 1000000 ;
56 Eigen::Matrix<double,Eigen::Dynamic,1> x, y ;
57 generate_table( x, y, N ) ;
58
59 // define a matrix od double (d) of size 2x2
60 // Matrix2d = Matrix<double,2,2>
61 Eigen::Matrix2d A ;
62 Eigen::Vector2d b, sol ;
63
64 // sum xi*xi
65 A(0,0) = x.dot(x) ;
66
67 // sum xi
68 A(1,0) = A(0,1) = x.sum() ;
69 A(1,1) = N ;
70
71 // sum xi*yi
72 b(0) = x.dot(y) ;
73
74 // sum yi
75 b(1) = y.sum() ;
76
77 // use Cholewsky decomposition
78 Eigen::LLT<Eigen::Matrix2d> lltOfA(A) ;
79 sol = lltOfA.solve(b) ;
80 cout << "Solution = " << sol << '\n' ;
81
82 // solve the problem using least squares
83 Eigen::Matrix<double,Eigen::Dynamic,Eigen::Dynamic> Mat(N,2) ;
84 Eigen::Matrix<double,Eigen::Dynamic,1> Rhs(N) ;
85 Mat.col(0) = x ;
86 Mat.col(1).fill(1) ; // fill column with one's
87 Rhs = y ;
88
89 cout << "Mat^T * Mat =\n" << Mat.transpose() * Mat << '\n' ;
90 cout << "A =\n" << A << '\n' ;
91
92 sol = Mat.colPivHouseholderQr().solve(Rhs) ;
93
94 cout << "Solution (LS) = " << sol << '\n' ;
95
96 }
97 catch ( exception const & error ) {
98 cerr << "Exception number " << error.what() << "\n" ;
99 }
100 catch ( ... ) {
101 cerr << "Exception found: Unknown error\n" ;
102 }
103 cout << "END OF PROGRAM\n" ;
104
105 return 0 ;
106}
File example20.cc
:open:
1/*
2 Some experiments using Eigen3
3
4 -I/usr/local/include/eigen3 to add search of headers for eigen3 library
5
6
7 Linear regression, bad way and good way.
8 Using more Eigen3
9
10*/
11
12#include <iostream>
13#include <iomanip>
14#include <vector>
15#include <cmath>
16#include <string> // C++ string
17#include <exception> // std::exception
18
19// include headers for the Dense linear algebra of Eigen library
20#include <Eigen/Dense>
21
22using namespace std ;
23//using namespace Eigen ;
24
25template <typename T>
26T
27func( T x )
28{ return 10*x-3 ; }
29
30template <typename T>
31void
32generate_table( Eigen::Matrix<T,Eigen::Dynamic,1> & x,
33 Eigen::Matrix<T,Eigen::Dynamic,1> & y,
34 int N ) {
35 x.resize(N) ;
36 y.resize(N) ;
37
38 for ( int i = 0 ; i < N ; ++i ) {
39 T xx = ((rand() % 10000)-5000) * 0.001 ;
40 T yy = func(xx) ;
41 yy += ((rand() % 10000)-5000) * 0.0001 ;
42 x(i) = xx ;
43 y(i) = yy ;
44 }
45
46}
47
48typedef double valueType ;
49
50int
51main() {
52
53 try {
54 int N = 1000000 ;
55 Eigen::Matrix<double,Eigen::Dynamic,1> x, y ;
56 generate_table( x, y, N ) ;
57
58 // define a matrix od double (d) of size 2x2
59 // Matrix2d = Matrix<double,2,2>
60 Eigen::Matrix2d A ;
61 Eigen::Vector2d b, sol ;
62
63 valueType tmp = x.sum() ;
64 A << x.dot(x), tmp,
65 tmp, N ;
66
67 b << x.dot(y), y.sum() ;
68
69 // use Cholewsky decomposition
70 Eigen::LLT<Eigen::Matrix2d> lltOfA(A) ;
71 sol = lltOfA.solve(b) ;
72 cout << "Solution = " << sol << '\n' ;
73
74 // solve the problem using least squares
75 Eigen::Matrix<double,Eigen::Dynamic,Eigen::Dynamic> Mat(N,2) ;
76 Mat.col(0) = x ;
77 Mat.col(1).fill(1) ; // fill column by ones!
78
79 sol = Mat.colPivHouseholderQr().solve(y) ;
80
81 cout << "Solution (LS) = " << sol << '\n' ;
82
83 }
84 catch ( exception const & error ) {
85 cerr << "Exception number " << error.what() << "\n" ;
86 }
87 catch ( ... ) {
88 cerr << "Exception found: Unknown error\n" ;
89 }
90 cout << "END OF PROGRAM\n" ;
91
92 return 0 ;
93}
All the files in one zip file