\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
NA 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 AS EXPOSED ON
http://eigen.tuxfamily.org/index.php?title=Benchmark
\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. So we have handwritten code for e.g. matrix product.

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}