\documentclass{slides}

\usepackage[utf8]{inputenc}
\usepackage{amssymb, amsmath, amsthm, comment}
\usepackage[totalwidth=520pt]{geometry}
\hyphenpenalty=5000
\tolerance=1000

\begin{document}
\begin{slide}
\begin{center}
Eigen: a practical approach to linear algebra\\
http://eigen.tuxfamily.org
\vskip 2 cm
Benoît Jacob, Dept. of Mathematics,\\
bjacob@math.toronto.edu
\vskip 2 cm
CS seminar, University of Toronto\\
23 january 2009
\vskip 2 cm
Eigen is co-developed with\\
Gaël Guennebaud (INRIA Bordeaux)\\
and a handful of occasional contributors
\end{center}
\end{slide}

\begin{slide}
Eigen:
\begin{itemize}
\item is a C++ library for linear algebra.
\item is versatile, all-in-one
\item has a great C++ API
\item has great performance
\item is lightweight
\item is LGPL licensed and cross-platform
\item is used in real world projects
\end{itemize}
\end{slide}

\begin{slide}
Eigen is {\bf versatile}

All 3 kinds of matrices/vectors:
\begin{itemize}
\item Dense, fixed-size
\item Dense, dynamic-size
\item Sparse (experimental)
\end{itemize}
Moreover these 3 kinds are fully {\bf integrated} with one another.

The fixed-size case is crucial:
\begin{itemize}
\item Very commonly used
\item Huge optimization opportunity
\end{itemize}
\end{slide}

\begin{slide}
Eigen is {\bf all-in-one}
\begin{itemize}
\item Basic matrix/vector functionality
\item Linear algebra algorithms (LU, QR, SVD, ...)
\item Geometry framework (projective, quaternions, ...)
\item Array manipulation
\end{itemize}

Fully {\bf self-contained} for {\bf dense} algorithms.

Allows optional {\bf backends} for {\bf sparse} algorithms.
\end{slide}

\begin{slide}
Eigen has a {\bf great C++ API}

We'll give many examples shortly.

Teaser:

{\tt\begin{verbatim}
matrix.row(i) += lambda * matrix.row(j);
\end{verbatim}}

Notice: in C++, a+=b means a=a+b.

Eigen works on expressions and decides\\{\bf automatically} when to use lazy evaluation.
\end{slide}

\begin{slide}
Eigen has {\bf great performance}
\begin{itemize}
\item Works at the level of expressions
\item Explicit vectorization: SSE2+ and AltiVec
\item Cache-friendly algorithms
\item Takes advantage of fixed dimensions
\end{itemize}

We'll show some benchmarks...
\end{slide}

\begin{slide}
Eigen is {\bf lightweight}

It is a {\bf pure template} library.

Compiled directly into the user code.

Only the code that's {\bf actually used}, is compiled.

Compilation times are still reasonable.

No binary library to link to.

Eigen itself is small: 16,000 LOC.
\end{slide}

\begin{slide}
Eigen is {\bf available}

License: LGPL 3+ (alternatively GPL 2+)

Closed-source software may use Eigen.

Supported operating systems:\\
Linux, BSD, MacOSX, Windows

Supported compilers:\\
GCC 3.3+, MSVC 2005+, ICC

For vectorization:\\
Instruction set: SSE2+, AltiVec\\
Compiler: GCC 4.2+, MSVC 2008+, ICC
\end{slide}

\begin{slide}
Eigen is used in real world apps:\\
\begin{itemize}
\item Computer graphics: MeshLab, VcgLib, Krita, libmv
\item Robotics companies: Yujin, Willow Garage
\item Desktop: KDE, KOffice (including Krita)
\item Chemistry: Avogadro, soon Open Babel
\end{itemize}
\end{slide}

\begin{slide}
Example program:

{\tt\begin{verbatim}
#include <Eigen/Core>
using namespace Eigen;
using namespace std;

int main()
{
  Matrix4f m = Matrix4f::Zero();
  m.diagonal().end(2) << 1, 2;
  cout << m << endl;
}
\end{verbatim}}

Matrix4f is a shortcut for:
{\tt\begin{verbatim}
Matrix<float,4,4>
\end{verbatim}}
\end{slide}

\begin{slide}
Matlab code:
{\tt\begin{verbatim}
m = eye(20);
m = rand(20);
\end{verbatim}}

Eigen/C++ equivalent:
{\tt\begin{verbatim}
MatrixXf m = MatrixXf::Identity(20,20);
MatrixXf m = MatrixXf::Random(20,20);
\end{verbatim}}

MatrixXf is a shortcut for:
{\tt\begin{verbatim}
Matrix<float,Dynamic,Dynamic>
\end{verbatim}}
\end{slide}

\begin{slide}
Matlab code:
{\tt\begin{verbatim}
v = rand(20, 1);
w = rand(20, 1); w = w / norm(w);
m = diag(rand(20, 1));
\end{verbatim}}

Eigen/C++ equivalent:
{\tt\begin{verbatim}
VectorXf v = VectorXf::Random(20);
VectorXf w = VectorXf::Random(20).normalized();
MatrixXf m = VectorXf::Random(20).asDiagonal();
\end{verbatim}}

VectorXf is a shortcut for:
{\tt\begin{verbatim}
Matrix<float,Dynamic,1>
\end{verbatim}}
\end{slide}

\begin{slide}
Algebraic expressions work as expected
{\tt\begin{verbatim}
m4 = m3 * (x1 * m1 + x2 * m2);
\end{verbatim}}
and produce optimized code.

No bad aliasing effects:
{\tt\begin{verbatim}
m = m * m;
\end{verbatim}}

Eigen is conservative with operator overloading.
\end{slide}

\begin{slide}
Matlab code:
{\tt\begin{verbatim}
m(i, j) = 0;
m = zeros(size(m));
m(i, :) = 0;
\end{verbatim}}

Eigen/C++ equivalent:
{\tt\begin{verbatim}
m(i,j) = 0;
m.setZero();
m.row(i).setZero();
\end{verbatim}}
\end{slide}

\begin{slide}
Matlab code:
{\tt\begin{verbatim}
m(i, :) = m(i, :) + x*m(j, :);
m(:, [i j]) = m(:, [j i]);
\end{verbatim}}

Eigen/C++ equivalent:
{\tt\begin{verbatim}
m.row(i) += x * m.row(j);
m.col(i).swap(m.col(j));
\end{verbatim}}
\end{slide}

\begin{slide}
Matlab code:
{\tt\begin{verbatim}
m(i:i+n-1, j:j+n-1) = eye(n);
m(i:i+n-1, j:j+n-1) = m2(i:i+n-1, j:j+n-1) * m3;
\end{verbatim}}

Eigen/C++ equivalent:
{\tt\begin{verbatim}
m.block(i,j,n,n).setIdentity();
m.block(i,j,n,n) = m2.block(i,j,n,n) * m3;
\end{verbatim}}
\end{slide}

\begin{slide}
Matlab code:
{\tt\begin{verbatim}
result = sum(m(:));
result = sum(m(:, i));
result = sum(m');
result = sum(m(i, :).^3);
\end{verbatim}}

Eigen/C++ equivalent:
{\tt\begin{verbatim}
result = m.sum();
result = m.col(i).sum();
result = m.rowwise().sum();
result = m.row(i).cwise().cube().sum();
\end{verbatim}}
\end{slide}

\begin{slide}
Matlab code:
{\tt\begin{verbatim}
m = [eye(n, p),    zeros(n, p);
     random(n, p), m2+m3];

v(i:i+3) = [a b c d];
\end{verbatim}}

Eigen/C++ equivalent:
{\tt\begin{verbatim}
m << MatrixXf::Identity(n,p), MatrixXf::Zero(n,p),
     MatrixXf::Random(n,p),   m2+m3;

v.segment(i,4) << a,b,c,d;
\end{verbatim}}
\end{slide}

\begin{slide}
\begin{center}
(SHOW BENCHMARKS)
\end{center}
\end{slide}

\begin{slide}
Eigen internals

With Eigen, an operation like
{\tt\begin{verbatim}
v3 = v1 + v2;
\end{verbatim}}
Is entirely performed in the "=".

The "+" by itself does nothing. It just returns a "sum expression" object.

The whole expression tree is encoded in the {\bf type} of this object.

Example: the expression 
{\tt\begin{verbatim}
m1 + m2.transpose()
\end{verbatim}}
gives an object of type
{\tt\begin{verbatim}
Sum< MatrixXf, Transpose<MatrixXf> >
\end{verbatim}}
\end{slide}

\begin{slide}
The expression types carry a lot of recursive {\bf metadata} about the expressions:

$\bullet$ Scalar type\\
$\bullet$ Dimensions at compile-time, or Dynamic\\
$\bullet$ Maximum dimensions at compile-time, or Dynamic\\
$\bullet$ Linearity: index-based access\\
$\bullet$ Cost of reading a coefficient\\
$\bullet$ Hints for/against lazy evaluation\\
$\bullet$ Vectorizability: packet access\\
$\bullet$ Data alignment\\
$\bullet$ Prefered storage order\\
$\bullet$ Shape (diagonal...)\\
\end{slide}

\begin{slide}
By default we {\bf lazy evaluate} expression.

However this is not always good.

Eigen automatically takes the decision to evaluate an intermediate step into a temporary.

Eigen then takes advantage of that to split the expression tree, limiting its depth.
\end{slide}

\begin{slide}
Lazy Evaluation:

$\bullet$ Can be {\bf necessary}.\\
Example: m1.block(r,c,nr,nc) = m2.\\
$\bullet$ Can be {\bf very good}.\\
Example: v4=v1+v2+v3.\\
$\bullet$ Can be {\bf bad}.\\
Example: m4=(m1+m2)*m3\\
$\bullet$ Can be {\bf terrible}.\\
Example: m4=(m1*m2)*m3\\
$\bullet$ Can be {\bf dangerous}.\\
Example: m = m*m

We conducted {\bf experiments}.

Lazy evaluation must be justified by a {\bf clear} gain.
\end{slide}

\begin{slide}
{\bf Vectorization}

Dynamic-size matrices: we allow ourselves some runtime logic (taking care of unaligned boundaries...)

Fixed-size matrices: no runtime overhead allowed

The expression metadata allows to choose the right vectorization strategy.

When cache friendliness matters, nothing replaces handwritten code.

Portable abstractions for SIMD instructions and packets.
\end{slide}

\begin{slide}
Vectorization Example.

{\tt\begin{verbatim}
#include<Eigen/Core>
using namespace Eigen;

void foo(Matrix2f& u,
         float a, const Matrix2f& v,
         float b, const Matrix2f& w)
{
  u = a*v + b*w - u;
}
\end{verbatim}}
\end{slide}

\begin{slide}
Comments:
\begin{itemize}
\item Neither rows not columns fit in 128bit packets
\item So necessary to understand that the expression is {\bf linearizable}
\item Notice how lazy evaluation allows to do only 1 traversal
\end{itemize}
\end{slide}

\begin{slide}
{\tt\begin{verbatim}
movl  8(%ebp), %edx
movss 20(%ebp), %xmm0
movl  24(%ebp), %eax
movaps  %xmm0, %xmm2
shufps  $0, %xmm2, %xmm2
movss 12(%ebp), %xmm0
movaps  %xmm2, %xmm1
mulps (%eax), %xmm1
shufps  $0, %xmm0, %xmm0
movl  16(%ebp), %eax
mulps (%eax), %xmm0
addps %xmm1, %xmm0
subps (%edx), %xmm0
movaps  %xmm0, (%edx)
\end{verbatim}}
\end{slide}

\begin{slide}
\begin{center}
{\Large Comments? Questions?}
\end{center}

http://eigen.tuxfamily.org

bjacob@math.toronto.edu

{\bf Core developers:} Ga\"el Guennebaud, Beno\^it Jacob

{\bf Contributors:} David Benjamin, Armin Berres, Daniel G\'omez, Konstantinos Margaritis, Christian Mayer,
Keir Mierle, Michael Olbrich, Kenneth Riddile, Alex Stapleton

\end{slide}

\end{document}
