Understanding Eigen: High-Level C++ Library for Linear Algebra

Eigen is a high-level C++ library offering a range of functionalities for linear algebra, matrix and vector operations, geometrical transformations, numerical solvers, and related algorithms. It provides efficient multidimensional array storage, fast math operations, and linear algebra capabilities. By leveraging Eigen, users can avoid cumbersome data splitting and easily manipulate data objects while retaining their native dimensionality. This article delves into Eigen's features, advantages, and practical usage through code examples and explanations.

• Eigen library
• C++
• linear algebra
• matrix operations
• numerical solvers

carson Follow

Uploaded on Jul 16, 2024 | 0 Views

Download Presentation

Please find below an Image/Link to download the presentation.

The content on the website is provided AS IS for your information and personal use only. It may not be sold, licensed, or shared on other websites without obtaining consent from the author. Download presentation by click this link. If you encounter any issues during the download, it is possible that the publisher has removed the file from their server.

Presentation Transcript

1. Eigen Ryan Honeyager honeyage@ucar.edu March 12, 2020

2. What is Eigen? http://eigen.tuxfamily.org/ Eigen is a high-level C++ library of template headers for linear algebra, matrix and vector operations, geometrical transformations, numerical solvers and related algorithms. We want to use Eigen for three purposes: Multidimensional array storage Very fast math Linear algebra 2

3. Why Do We Want It? We want our data objects to retain their native dimensionality. It is much easier to read and perform coefficientwise math without looping over index variables. We want to avoid situations where we split up data (e.g. Tb_ch1, Tb_ch2, ). Eigen provides classes for manipulating data that are very efficient and easy to use. Eigen is already a dependency of OOPS. We want to extend the IODA ObsSpace / ObsData / ObsDataVector interfaces to natively store and manipulate Eigen objects. 3

4. A Very Simple Program CMakeLists.txt #include <iostream> #include <Eigen/Dense> using Eigen::ArrayXXd; cmake_minimum_required(VERSION 3.2) project(ex VERSION 0.0.1 LANGUAGES C CXX) find_package(Eigen3 REQUIRED) int main(int, char**) { ArrayXXd a(2,2); a(0,0) = 3; a(1,0) = 2.5; a(0,1) = -1; a(1,1) = a(1,0) + a(0,1); std::cout << a << std::endl; } Create a 2x2 array. add_executable(ex1 ex1.cpp) set_property(TARGET ex1 PROPERTY CXX_STANDARD 14) target_include_directories(ex1 SYSTEM PUBLIC ${EIGEN3_INCLUDE_DIR}) Set the array coefficients. Print the array. 3 -1 2.5 1.5 4

5. What is an Eigen::Array? Dynamic allocation allows objects to be resized. If you have something reasonably small, you can try to allocate on the stack instead. Other types: char float int Eigen Arrays are templates. Storage Type Rows and columns at compile time. Generally should be Eigen::Dynamic. Eigen::ArrayXXd Eigen::Array< double, Eigen::Dynamic, Eigen::Dynamic, Eigen::ColMajor> 1 4 7 2 5 8 3 6 9 Store in Column-Major Form Row-major: Column-major: 1, 2, 3, 4, 5, 6, 7, 8, 9 1, 4, 7, 2, 5, 8, 3, 6, 9 5

6. Math with Eigen Eigen::Arrays can perform coefficientwise math. So, you can add, subtract, multiply, divide. ArrayXXd a(2,2), b(2,2); a << 1, 2, 3, 4; b << 2, 3, 4, 5; ? =1 ? =2 ? =3 ? =4 2 4 3 5 5 9 5 7 3 4 auto c = a + b; auto d = ((a*2) b) + 4; 7 6 Eigen::Matrix is another Eigen class that stores data like an array, but implements matrix math (matrix multiplication, outer products, determinants, etc.). 6

7. How Eigen Works Behind The Scenes Eigen is implemented using the expression templates metaprogramming technique, meaning it builds expression trees at compile time and generates custom code to evaluate these. Using expression templates and a cost model of floating point operations, the library performs its own loop unrolling and vectorization. http://eigen.tuxfamily.org/dox/TopicInsideEigenExample.html 1 for(int i = 0; i < size; i++) tmp[i] = v[i] + w[i]; for(int i = 0; i < size; i++) u[i] = tmp[i]; #include<Eigen/Core> int main() { int size = 50; // VectorXf is a vector of floats, with dynamic size. Eigen::VectorXf u(size), v(size), w(size); u = v + w; } 3 4 for(int i = 0; i < 4*(size/4); i+=4) u.packet(i) = v.packet(i) + w.packet(i); for(int i = 4*(size/4); i < size; i++) u[i] = v[i] + w[i]; Explanations: 1. Adding two vectors to get a third. 2. This gets implemented elementwise. 3. In a na ve implementation of this class, operator+ returns a temporary object, which we then have to copy to u. Eigen uses lazy evaluation to simplify algebraic logic at compile time, and takes advantage of return value optimization. 4. Eigen knows that it can vectorize the loop into groups of 4 to use SSE2 and AltiVec instructions for speed. 2 for(int i = 0; i < size; i++) u[i] = v[i] + w[i]; 7

8. Manipulating Array Objects A thorough tutorial is provided in the Eigen documentation: http://eigen.tuxfamily.org/dox/group__TutorialMatrixClass.html ArrayXXd a(2,2); a << 1, 2, 3, 4; Another way to set initial values Resizing a.resize(6,8); Resize a to a 6x8 array. Block Operations a.block(2,0,4,8); Start row, start col, # rows, # cols Converting data types auto b = a.cast<float>(); Setting to a constant value a.setConstant(0); Convert an array of doubles to an array of floats. Sets all of a to a constant zero. a.block(2,0,4,8).setConstant(3); Sets the block to 3. 8

9. Reductions with Eigen ArrayXXd a(2,2); a << 1, 2, 3, 4; cout << "a.sum(): " << a.sum() << endl; cout << "a.prod(): " << a.prod() << endl; cout << "a.mean(): " << a.mean() << endl; cout << "a.minCoeff(): " << a.minCoeff() << endl; cout << "a.maxCoeff(): " << a.maxCoeff() << endl; cout << "a.trace(): " << a.trace() << endl; a.sum(): 10 a.prod(): 24 a.mean(): 2.5 a.minCoeff(): 1 a.maxCoeff(): 4 a.trace(): 5 9

10. Eigen::Map Getting Large Data into Eigen Occasionally you may have a pre-defined array of numbers that you want to use within Eigen as a vector or matrix. While one option is to make a copy of the data, most commonly you probably want to re-use this memory as an Eigen type. An Eigen::Map is a view . It treats data as an Eigen object without making a copy. Map<Matrix<typename Scalar, int RowsAtCompileTime, int ColsAtCompileTime> > int in_data[6] = {0, 1, 2, 3, 4, 5}; int rows = 2, cols = 3; Map<ArrayXXi> mapped_in_data(in_data,rows,columns); ArrayXXi some_math_result = (mapped_in_data * 4) + 6; 10

11. Getting Data Out of Eigen You might want to convert Eigen objects to 1-D arrays. Use the .data() method. Eigen::ArrayXXf some_floats(3,4); some_floats << 0.5f, 1.5f, 2.5f, 3.5f, ; std::vector<float> data_vector(some_floats.data(), some_floats.size()); Note: the .data() method is a raw data view. Row-wise and column-wise storage use different element ordering. 1 2 3 4 5 6 7 8 9 Row-major: Column-major: 1, 2, 3, 4, 5, 6, 7, 8, 9 1, 4, 7, 2, 5, 8, 3, 6, 9 11

12. What About Multidimensional Arrays? Eigen::Arrays are good for one and two-dimensional data. We still want to avoid using constructs like vector<vector<vector<vector<double>>>>. For higher dimensionalities, we can use Eigen::Tensor. Eigen::Tensor<data_type, rank>(size_0, size_1, ) // Map a tensor of ints on top of stack-allocated storage. int storage[128]; // 2 x 4 x 2 x 8 = 128 TensorMap<Tensor<int, 4>> t_4d(storage, 2, 4, 2, 8); Eigen::Tensors are still in development, but are stable as data containers. #include <unsupported/Tensor> 12