Numerical Linear Algebra and Factorizations¶

Contents¶

• Numerical Linear Algebra and Factorizations
  • Overview
  • Factorizations
  • Continuous-Time Markov Chains (CTMCs)
  • Banded Matrices
  • Implementation Details and Performance
  • Exercises

You cannot learn too much linear algebra. – Benedict Gross

Overview¶

In this lecture, we examine the structure of matrices and linear operators (e.g., dense, sparse, symmetric, tridiagonal, banded) and discuss how the structure can be exploited to radically increase the performance of solving large problems.

We build on applications discussed in previous lectures: linear algebra, orthogonal projections, and Markov chains.

The methods in this section are called direct methods, and they are qualitatively similar to performing Gaussian elimination to factor matrices and solve systems of equations. In iterative methods and sparsity we examine a different approach, using iterative algorithms, where we can think of more general linear operators.

The list of specialized packages for these tasks is enormous and growing, but some of the important organizations to look at are JuliaMatrices , JuliaSparse, and JuliaMath

NOTE: As this section uses advanced Julia techniques, you may wish to review multiple-dispatch and generic programming in introduction to types, and consider further study on generic programming.

The theme of this lecture, and numerical linear algebra in general, comes down to three principles:

1. Identify structure (e.g., symmetric, sparse, diagonal) matrices in order to use specialized algorithms.
2. Do not lose structure by applying the wrong numerical linear algebra operations at the wrong times (e.g., sparse matrix becoming dense)
3. Understand the computational complexity of each algorithm, given the structure of the inputs.

Setup¶

Computational Complexity¶

Ask yourself whether the following is a computationally expensive operation as the matrix size increases

• Multiplying two matrices?

  • Answer: It depends. Multiplying two diagonal matrices is trivial.

• Solving a linear system of equations?

  • Answer: It depends. If the matrix is the identity, the solution is the vector itself.

• Finding the eigenvalues of a matrix?

  • Answer: It depends. The eigenvalues of a triangular matrix are the diagonal elements.

As the goal of this section is to move toward numerical methods with large systems, we need to understand how well algorithms scale with the size of matrices, vectors, etc. This is known as computational complexity. As we saw in the answer to the questions above, the algorithm - and hence the computational complexity - changes based on matrix structure.

While this notion of complexity can work at various levels, such as the number of significant digits for basic mathematical operations, the amount of memory and storage required, or the amount of time, we will typically focus on the time complexity.

For time complexity, the size $ N $ is usually the dimensionality of the problem, although occasionally the key will be the number of non-zeros in the matrix or the width of bands. For our applications, time complexity is best thought of as the number of floating point operations (e.g., addition, multiplication) required.

Notation¶

Complexity of algorithms is typically written in Big O notation, which provides bounds on the scaling of the computational complexity with respect to the size of the inputs.

Formally, if the number of operations required for a problem size $ N $ is $ f(N) $, we can write this as $ f(N) = O(g(N)) $ for some $ g(N) $ - typically a polynomial.

The interpretation is that there exist some constants $ M $ and $ N_0 $ such that

$$ f(N) \leq M g(N), \text{ for } N > N_0 $$

For example, the complexity of finding an LU Decomposition of a dense matrix is $ O(N^3) $, which should be read as there being a constant where eventually the number of floating point operations required to decompose a matrix of size $ N\times N $ grows cubically.

Keep in mind that these are asymptotic results intended for understanding the scaling of the problem, and the constant can matter for a given fixed size.

For example, the number of operations required for an LU decomposition of a dense $ N \times N $ matrix is $ f(N) = \frac{2}{3} N^3 $, ignoring the $ N^2 $ and lower terms. Other methods of solving a linear system may have different constants of proportionality, even if they have the same scaling, $ O(N^3) $.

Rules of Computational Complexity¶

You will sometimes need to think through how combining algorithms changes complexity. For example, if you use

1. an $ O(N^3) $ operation $ P $ times, then it simply changes the constant. The complexity remains $ O(N^3) $.
2. one $ O(N^3) $ operation and one $ O(N^2) $ operation, then you take the max. The complexity remains $ O(N^3) $.
3. a repetition of an $ O(N) $ operation that itself uses an $ O(N) $ operation, you take the product. The complexity becomes $ O(N^2) $.

With this, we have an important word of caution: Dense-matrix multiplication is an expensive operation for unstructured matrices. The naive version is $ O(N^3) $ while the fastest-known algorithms (e.g., Coppersmith-Winograd) are roughly $ O(N^{2.37}) $. In practice, it is reasonable to crudely approximate with $ O(N^3) $ when doing an analysis, in part since the higher constant factors of the better scaling algorithms dominate the better complexity until matrices become very large.

Of course, modern libraries use highly tuned and numerically stable algorithms to multiply matrices and exploit the computer architecture, memory cache, etc., but this simply lowers the constant of proportionality and they remain roughly approximated by $ O(N^3) $.

A consequence is that, since many algorithms require matrix-matrix multiplication, it is often not possible to go below that order without further matrix structure.

That is, changing the constant of proportionality for a given size can help, but in order to achieve better scaling you need to identify matrix structure (e.g., tridiagonal, sparse) and ensure that your operations do not lose it.

Losing Structure¶

As a first example of a structured matrix, consider a sparse array.

This increase from less than 50 to 100 percent dense demonstrates that significant sparsity can be lost when computing an inverse.

The results can be even more extreme. Consider a tridiagonal matrix of size $ N \times N $ that might come out of a Markov chain or a discretization of a diffusion process,

The number of non-zeros here is approximately $ 3 N $, linear, which scales well for huge matrices into the millions or billions

But consider the inverse

Now, the matrix is fully dense and has $ N^2 $ non-zeros.

This also applies to the $ A' A $ operation when forming the normal equations of linear least squares.

We see that a 30 percent dense matrix becomes almost full dense after the product is taken.

Sparsity/Structure is not just for storage: Matrix size can sometimes become important (e.g., a 1 million by 1 million tridiagonal matrix needs to store 3 million numbers (i.e., about 6MB of memory), where a dense one requires 1 trillion (i.e., about 1TB of memory)).

But, as we will see, the main purpose of considering sparsity and matrix structure is that it enables specialized algorithms, which typically have a lower computational order than unstructured dense, or even unstructured sparse, operations.

First, create a convenient function for benchmarking linear solvers

Then, take away structure to see the impact on performance,

This example shows what is at stake: using a structured tridiagonal matrix may be 10-20 times faster than using a sparse matrix, which is 100 times faster than using a dense matrix.

In fact, the difference becomes more extreme as the matrices grow. Solving a tridiagonal system is $ O(N) $, while that of a dense matrix without any structure is $ O(N^3) $. The complexity of a sparse solution is more complicated, and scales in part by the nnz(N), i.e., the number of nonzeros.

Matrix Multiplication¶

While we write matrix multiplications in our algebra with abundance, in practice the computational operation scales very poorly without any matrix structure.

Matrix multiplication is so important to modern computers that the constant of scaling is small using proper packages, but the order is still roughly $ O(N^3) $ in practice (although smaller in theory, as discussed above).

Sparse matrix multiplication, on the other hand, is $ O(N M_A M_B) $ where $ M_A $ is the number of nonzeros per row of $ A $ and $ M_B $ is the number of non-zeros per column of $ B $.

By the rules of computational order, that means any algorithm requiring a matrix multiplication of dense matrices requires at least $ O(N^3) $ operation.

The other important question is what is the structure of the resulting matrix. For example, multiplying an upper triangular matrix by a lower triangular matrix

But the product is fully dense (e.g., think of a Cholesky multiplied by itself to produce a covariance matrix)

On the other hand, a tridiagonal matrix times a diagonal matrix is still tridiagonal - and can use specialized $ O(N) $ algorithms.

Factorizations¶

When you tell a numerical analyst you are solving a linear system using direct methods, their first question is “which factorization?”.

Just as you can factor a number (e.g., $ 6 = 3 \times 2 $) you can factor a matrix as the product of other, more convenient matrices (e.g., $ A = L U $ or $ A = Q R $, where $ L, U, Q, $ and $ R $ have properties such as being triangular, orthogonal, etc.).

Inverting Matrices¶

On paper, since the Invertible Matrix Theorem tells us that a unique solution is equivalent to $ A $ being invertible, we often write the solution to $ A x = b $ as

$$ x = A^{-1} b $$

What if we do not (directly) use a factorization?

Take a simple linear system of a dense matrix,

On paper, we try to solve the system $ A x = b $ by inverting the matrix,

As we will see throughout, inverting matrices should be used for theory, not for code. The classic advice that you should never invert a matrix may be slightly exaggerated, but is generally good advice.

Solving a system by inverting a matrix is always a little slower, is potentially less accurate, and will sometimes lose crucial sparsity compared to using factorizations. Moreover, the methods used by libraries to invert matrices are frequently the same factorizations used for computing a system of equations.

Even if you need to solve a system with the same matrix multiple times, you are better off factoring the matrix and using the solver rather than calculating an inverse.

Triangular Matrices and Back/Forward Substitution¶

Some matrices are already in a convenient form and require no further factoring.

For example, consider solving a system with an UpperTriangular matrix,

This system is especially easy to solve using back substitution. In particular, $ x_3 = b_3 / U_{33}, x_2 = (b_2 - x_3 U_{23})/U_{22} $, etc.

A LowerTriangular matrix has similar properties and can be solved with forward substitution.

The computational order of back substitution and forward substitution is $ O(N^2) $ for dense matrices. Those fast algorithms are a key reason that factorizations target triangular structures.

LU Decomposition¶

The $ LU $ decomposition finds a lower triangular matrix $ L $ and an upper triangular matrix $ U $ such that $ L U = A $.

For a general dense matrix without any other structure (i.e., not known to be symmetric, tridiagonal, etc.) this is the standard approach to solve a system and exploit the speed of back and forward substitution using the factorization.

The computational order of LU decomposition itself for a dense matrix is $ O(N^3) $ - the same as Gaussian elimination - but it tends to have a better constant term than others (e.g., half the number of operations of the QR decomposition). For structured or sparse matrices, that order drops.

We can see which algorithm Julia will use for the \ operator by looking at the factorize function for a given matrix.

In this case, it provides an $ L $ and $ U $ factorization (with pivoting ).

With the factorization complete, we can solve different b right hand sides.

In practice, the decomposition also includes a $ P $ which is a permutation matrix such that $ P A = L U $.

We can also directly calculate an LU decomposition with lu but without the pivoting,

And we can verify the decomposition

To see roughly how the solver works, note that we can write the problem $ A x = b $ as $ L U x = b $. Let $ U x = y $, which breaks the problem into two sub-problems.

$$ \begin{aligned} L y &= b\\ U x &= y \end{aligned} $$

As we saw above, this is the solution to two triangular systems, which can be efficiently done with back or forward substitution in $ O(N^2) $ operations.

To demonstrate this, first using

The LU decomposition also has specialized algorithms for structured matrices, such as a Tridiagonal

This factorization is the key to the performance of the A \ b in this case. For Tridiagonal matrices, the LU decomposition is $ O(N^2) $.

Finally, just as a dense matrix without any structure uses an LU decomposition to solve a system, so will the sparse solvers

With sparsity, the computational order is related to the number of non-zeros rather than the size of the matrix itself.

Cholesky Decomposition¶

For real, symmetric, positive semi-definite matrices, a Cholesky decomposition is a specialized example of an LU decomposition where $ L = U' $.

The Cholesky is directly useful on its own (e.g., Classical Control with Linear Algebra), but it is also an efficient factorization to use in solving symmetric positive semi-definite systems.

As always, symmetry allows specialized algorithms.

Here, the $ A $ decomposition is Bunch-Kaufman rather than Cholesky, because Julia doesn’t know that the matrix is positive semi-definite. We can manually factorize with a Cholesky,

Benchmarking,

QR Decomposition¶

Previously, we learned about applications of the QR decomposition to solving the linear least squares.

While in principle the solution to the least-squares problem

$$ \min_x \| Ax -b \|^2 $$

is $ x = (A'A)^{-1}A'b $, in practice note that $ A'A $ becomes dense and calculating the inverse is rarely a good idea.

The QR decomposition is a decomposition $ A = Q R $ where $ Q $ is an orthogonal matrix (i.e., $ Q'Q = Q Q' = I $) and $ R $ is an upper triangular matrix.

Given the previous derivation, we showed that we can write the least-squares problem as the solution to

$$ R x = Q' b $$

where, as discussed above, the upper-triangular structure of $ R $ can be solved easily with back substitution.

The \ operator solves the linear least-squares problem whenever the given A is rectangular

To manually use the QR decomposition in solving linear least squares:

This stacks the R with zeros, but the more specialized algorithm would not multiply directly in that way.

In some cases, if an LU is not available for a particular matrix structure, the QR factorization can also be used to solve systems of equations (i.e., not just LLS). This tends to be about 2 times slower than the LU but is of the same computational order.

Deriving the approach, where we can now use the inverse since the system is square and we assumed $ A $ was non-singular,

$$ \begin{aligned} A x &= b\\ Q R x &= b\\ Q^{-1} Q R x &= Q^{-1} b\\ R x &= Q' b \end{aligned} $$

where the last step uses the fact that $ Q^{-1} = Q' $ for an orthogonal matrix.

Given the decomposition, the solution for dense matrices is of computational order $ O(N^2) $. To see this, look at the order of each operation.

• Since $ R $ is an upper-triangular matrix, it can be solved quickly through back substitution with computational order $ O(N^2) $
• A transpose operation is of order $ O(N^2) $
• A matrix-vector product is also $ O(N^2) $

In all cases, the order would drop depending on the sparsity pattern of the matrix (and corresponding decomposition). A key benefit of a QR decomposition is that it tends to maintain sparsity.

Without implementing the full process, you can form a QR factorization with qr and then use it to solve a system

Spectral Decomposition¶

A spectral decomposition, also known as an eigendecomposition, finds all of the eigenvectors and eigenvalues to decompose a square matrix A such that

$$ A = Q \Lambda Q^{-1} $$

where $ Q $ is a matrix made of the eigenvectors of $ A $ as columns, and $ \Lambda $ is a diagonal matrix of the eigenvalues. Only square, diagonalizable matrices have an eigendecomposition (where a matrix is not diagonalizable if it does not have a full set of linearly independent eigenvectors).

In Julia, whenever you ask for a full set of eigenvectors and eigenvalues, it decomposes using an algorithm appropriate for the matrix type. For example, symmetric, Hermitian, and tridiagonal matrices have specialized algorithms.

To see this,

Keep in mind that a real matrix may have complex eigenvalues and eigenvectors, so if you attempt to check Q * Λ * inv(Q) - A - even for a positive-definite matrix - it may not be a real number due to numerical inaccuracy.

Continuous-Time Markov Chains (CTMCs)¶

In the previous lecture on discrete-time Markov chains, we saw that the transition probability between state $ x $ and state $ y $ was summarized by the matrix $ P(x, y) := \mathbb P \{ X_{t+1} = y \,|\, X_t = x \} $.

As a brief introduction to continuous time processes, consider the same state space as in the discrete case: $ S $ is a finite set with $ n $ elements $ \{x_1, \ldots, x_n\} $.

A Markov chain $ \{X_t\} $ on $ S $ is a sequence of random variables on $ S $ that have the Markov property.

In continuous time, the Markov Property is more complicated, but intuitively is the same as the discrete-time case.

That is, knowing the current state is enough to know probabilities for future states. Or, for realizations $ x(\tau)\in S, \tau \leq t $,

$$ \mathbb P \{ X(t+s) = y \,|\, X(t) = x, X(\tau) = x(\tau) \text{ for } 0 \leq \tau \leq t \} = \mathbb P \{ X(t+s) = y \,|\, X(t) = x\} $$

Heuristically, consider a time period $ t $ and a small step forward, $ \Delta $. Then the probability to transition from state $ i $ to state $ j $ is

$$ \mathbb P \{ X(t + \Delta) = j \,|\, X(t) \} = \begin{cases} q_{ij} \Delta + o(\Delta) & i \neq j\\ 1 + q_{ii} \Delta + o(\Delta) & i = j \end{cases} $$

where the $ q_{ij} $ are “intensity” parameters governing the transition rate, and $ o(\Delta) $ is little-o notation. That is, $ \lim_{\Delta\to 0} o(\Delta)/\Delta = 0 $.

Just as in the discrete case, we can summarize these parameters by an $ N \times N $ matrix, $ Q \in R^{N\times N} $.

Recall that in the discrete case every element is weakly positive and every row must sum to one. With continuous time, however, the rows of $ Q $ sum to zero, where the diagonal contains the negative value of jumping out of the current state. That is,

• $ q_{ij} \geq 0 $ for $ i \neq j $
• $ q_{ii} \leq 0 $
• $ \sum_{j} q_{ij} = 0 $

The $ Q $ matrix is called the intensity matrix, or the infinitesimal generator of the Markov chain. For example,

$$ Q = \begin{bmatrix} -0.1 & 0.1 & 0 & 0 & 0 & 0\\ 0.1 &-0.2 & 0.1 & 0 & 0 & 0\\ 0 & 0.1 & -0.2 & 0.1 & 0 & 0\\ 0 & 0 & 0.1 & -0.2 & 0.1 & 0\\ 0 & 0 & 0 & 0.1 & -0.2 & 0.1\\ 0 & 0 & 0 & 0 & 0.1 & -0.1\\ \end{bmatrix} $$

In the above example, transitions occur only between adjacent states with the same intensity (except for a ``bouncing back’’ of the bottom and top states).

Implementing the $ Q $ using its tridiagonal structure

Here we can use Tridiagonal to exploit the structure of the problem.

Consider a simple payoff vector $ r $ associated with each state, and a discount rate $ ρ $. Then we can solve for the expected present discounted value in a way similar to the discrete-time case.

$$ \rho v = r + Q v $$

or rearranging slightly, solving the linear system

$$ (\rho I - Q) v = r $$

For our example, exploiting the tridiagonal structure,

Note that this $ A $ matrix is maintaining the tridiagonal structure of the problem, which leads to an efficient solution to the linear problem.

The $ Q $ is also used to calculate the evolution of the Markov chain, in direct analogy to the $ ψ_{t+k} = ψ_t P^k $ evolution with the transition matrix $ P $ of the discrete case.

In the continuous case, this becomes the system of linear differential equations

$$ \dot{ψ}(t) = Q(t)^T ψ(t) $$

given the initial condition $ \psi(0) $ and where the $ Q(t) $ intensity matrix is allowed to vary with time. In the simplest case of a constant $ Q $ matrix, this is a simple constant-coefficient system of linear ODEs with coefficients $ Q^T $.

If a stationary equilibrium exists, note that $ \dot{ψ}(t) = 0 $, and the stationary solution $ ψ^{*} $ needs to satisfy

$$ 0 = Q^T ψ^{*} $$

Notice that this is of the form $ 0 ψ^{*} = Q^T ψ^{*} $ and hence is equivalent to finding the eigenvector associated with the $ \lambda = 0 $ eigenvalue of $ Q^T $.

With our example, we can calculate all of the eigenvalues and eigenvectors

Indeed, there is a $ \lambda = 0 $ eigenvalue, which is associated with the last column in the eigenvector. To turn that into a probability, we need to normalize it.

Multiple Dimensions¶

A frequent case in discretized models is dealing with Markov chains with multiple “spatial” dimensions (e.g., wealth and income).

After discretizing a process to create a Markov chain, you can always take the Cartesian product of the set of states in order to enumerate it as a single state variable.

To see this, consider states $ i $ and $ j $ governed by infinitesimal generators $ Q $ and $ A $.

This provides the combined Markov chain for the $ (i,j) $ process. To see the sparsity pattern,

To calculate a simple dynamic valuation, consider whether the payoff of being in state $ (i,j) $ is $ r_{ij} = i + 2j $

Solving the equation $ \rho v = r + L v $

The reshape helps to rearrange it back to being two-dimensional.

To find the stationary distribution, we calculate the eigenvalue and choose the eigenvector associated with $ \lambda=0 $ . In this case, we can verify that it is the last one.

Reshape this to be two-dimensional if it is helpful for visualization.

Irreducibility¶

As with the discrete-time Markov chains, a key question is whether CTMCs are reducible, i.e., whether states communicate. The problem is isomorphic to determining whether the directed graph of the Markov chain is strongly connected.

We can verify that it is possible to move between every pair of states in a finite number of steps with

Alternatively, as an example of a reducible Markov chain where states $ 1 $ and $ 2 $ cannot jump to state $ 3 $.

Banded Matrices¶

A tridiagonal matrix has 3 non-zero diagonals: the main diagonal, the first sub-diagonal (i.e., below the main diagonal), and also the first super-diagonal (i.e., above the main diagonal).

This is a special case of a more general type called a banded matrix, where the number of sub- and super-diagonals can be greater than 1. The total width of main-, sub-, and super-diagonals is called the bandwidth. For example, a tridiagonal matrix has a bandwidth of 3.

An $ N \times N $ banded matrix with bandwidth $ P $ has about $ N P $ nonzeros in its sparsity pattern.

These can be created directly as a dense matrix with diagm. For example, with a bandwidth of three and a zero diagonal,

Or as a sparse matrix,

Or directly using BandedMatrices.jl

There is also a convenience function for generating random banded matrices

And, of course, specialized algorithms will be used to exploit the structure when solving linear systems. In particular, the complexity is related to the $ O(N P_L P_U) $ for upper and lower bandwidths $ P $

The factorization algorithm uses a specialized LU decomposition for banded matrices.

Implementation Details and Performance¶

Recall the famous quote from Knuth: “97% of the time, premature optimization is the root of all evil. Yet we should not pass up our opportunities in that critical 3%.” The most common example of premature optimization is trying to use your own mental model of a compiler while writing your code, worried about the efficiency of code, and (usually incorrectly) second-guessing the compiler.

Concretely, the lessons in this section are:

1. Don’t worry about optimizing your code unless you need to. Code clarity is your first-order concern.
2. If you use other people’s packages, they can worry about performance and you don’t need to.
3. If you absolutely need that “critical 3%,” your intuition about performance is usually wrong on modern CPUs and GPUs, so let the compiler do its job.
4. Benchmarking (e.g., @btime) and profiling are the tools to figure out performance bottlenecks. If 99% of computing time is spent in one small function, then there is no point in optimizing anything else.
5. If you benchmark to show that a particular part of the code is an issue, and you can’t find another library that does a better job, then you can worry about performance.

You will rarely get to step 3, let alone step 5.

However, there is also a corollary: “don’t pessimize prematurely.” That is, don’t make choices that lead to poor performance without any tradeoff in improved code clarity. For example, writing your own algorithms when a high-performance algorithm exists in a package or Julia itself, or lazily making a matrix dense and carelessly dropping its structure.

Implementation Difficulty¶

Numerical analysts sometimes refer to the lowest level of code for basic operations (e.g., a dot product, matrix-matrix product, convolutions) as kernels.

That sort of code is difficult to write, and performance depends on the characteristics of the underlying hardware, such as the instruction set available on the particular CPU, the size of the CPU cache, and the layout of arrays in memory.

Typically, these operations are written in a BLAS library, organized into different levels. The levels roughly correspond to the computational order of the operations: BLAS Level 1 are $ O(N) $ operations such as linear products, Level 2 are $ O(N^2) $ operations such as matrix-vector products, and Level 3 are roughly $ O(N^3) $, such as general matrix-matrix products.

An example of a BLAS library is OpenBLAS, which is used by default in Julia, or the Intel MKL, which is used in Matlab (and in Julia if the MKL.jl package is installed).

On top of BLAS are LAPACK operations, which are higher-level kernels, such as matrix factorizations and eigenvalue algorithms, and are often in the same libraries (e.g., MKL has both BLAS and LAPACK functionality).

The details of these packages are not especially relevant, but if you are talking about performance, people will inevitably start discussing these different packages and kernels. There are a few important things to keep in mind:

1. Leave writing kernels to the experts. Even simple-sounding algorithms can be very complicated to implement with high performance.
2. Your intuition about performance of code is probably going to be wrong. If you use high quality libraries rather than writing your own kernels, you don’t need to use your intuition.
3. Don’t get distracted by the jargon or acronyms above if you are reading about performance.

Row- and Column-Major Ordering¶

There is a practical performance issue which may influence your code. Since memory in a CPU is linear, dense matrices need to be stored by either stacking columns (called column-major order) or rows.

The reason this matters is that compilers can generate better performance if they work in contiguous chunks of memory, and this becomes especially important with large matrices due to the interaction with the CPU cache. Choosing the wrong order when there is no benefit in code clarity is an example of premature pessimization. The performance difference can be orders of magnitude in some cases, and nothing in others.

One option is to use the functions that let the compiler choose the most efficient way to traverse memory. If you need to choose the looping order yourself, then you might want to experiment with going through columns first and going through rows first. Other times, let Julia decide, i.e., enumerate and eachindex will choose the right approach.

Julia, Fortran, and Matlab all use column-major order, while C/C++ and Python use row-major order. This means that if you find an algorithm written for C/C++/Python, you will sometimes need to make small changes if performance is an issue.

Digression on Allocations and In-place Operations¶

While we have usually not considered optimizing code for performance (and have focused on the choice of algorithms instead), when matrices and vectors become large we need to be more careful.

The most important thing to avoid are excess allocations, which usually occur due to the use of temporary vectors and matrices when they are not necessary. Sometimes those extra temporary values can cause enormous degradations in performance.

However, caution is suggested since excess allocations are never relevant for scalar values, and allocations frequently create faster code for smaller matrices/vectors since it can lead to better cache locality.

To see this, a convenient tool is the benchmarking

The ! on the f! is an informal way to say that the function is mutating, and the first argument (C here) is by convention the modified variable.

In the f! function, notice that the D is a temporary variable which is created, and then modified afterwards. But notice that since C is modified directly, there is no need to create the temporary D matrix.

This is an example of where an in-place version of the matrix multiplication can help avoid the allocation.

Note that in the output of the benchmarking, the f2! is non-allocating and is using the pre-allocated C variable directly.

Another example of this is solutions to linear equations, where for large solutions you may pre-allocate and reuse the solution vector.

However, if you benchmark carefully, you will see that this is sometimes slower. Avoiding allocations is not always a good idea - and worrying about it prior to benchmarking is premature optimization.

There are a variety of other non-allocating versions of functions. For example,

Finally, a common source of unnecessary allocations is when taking slices or portions of matrices. For example, the following allocates a new matrix B and copies the values.

To see that these are different matrices, note that

Instead of allocating a new matrix, you can take a view of a matrix, which provides an appropriate AbstractArray type that doesn’t allocate new memory with the @view matrix.

But again, you will often find that doing @view leads to slower code. Benchmark instead, and generally rely on it for large matrices and for contiguous chunks of memory (e.g., columns rather than rows).

Exercises¶

Exercise 1¶

This exercise is for practice on writing low-level routines (i.e., “kernels”), and to hopefully convince you to leave low-level code to the experts.

The formula for matrix multiplication is deceptively simple. For example, with the product of square matrices $ C = A B $ of size $ N \times N $, the $ i,j $ element of $ C $ is

$$ C_{ij} = \sum_{k=1}^N A_{ik} B_{kj} $$

Alternatively, you can take a row $ A_{i,:} $ and column $ B_{:, j} $ and use an inner product

$$ C_{ij} = A_{i,:} \cdot B_{:,j} $$

Note that the inner product in a discrete space is simply a sum, and has the same complexity as the sum (i.e., $ O(N) $ operations).

For a dense matrix without any structure and using a naive multiplication algorithm, this also makes it clear why the complexity is $ O(N^3) $: You need to evaluate it for $ N^2 $ elements in the matrix and do an $ O(N) $ operation each time.

For this exercise, implement matrix multiplication yourself and compare performance in a few permutations.

1. Use the built-in function in Julia (i.e., C = A * B, or, for a better comparison, the in-place version mul!(C, A, B), which works with pre-allocated data).
2. Loop over each $ C_{ij} $ by the row first (i.e., the i index) and use a for loop for the inner product.
3. Loop over each $ C_{ij} $ by the column first (i.e., the j index) and use a for loop for the inner product.
4. Do the same but use the dot product instead of the sum.
5. Choose your best implementation of these, and then for matrices of a few different sizes (N=10, N=1000, etc.), and compare the ratio of performance of your best implementation to the built-in BLAS library.

A few more hints:

• You can just use random matrices (e.g., A = rand(N, N)).
• For all of them, pre-allocate the $ C $ matrix beforehand with C = similar(A) or something equivalent.
• To compare performance, put your code in a function and use the @btime macro to time it.

Exercise 2a¶

Here we will calculate the evolution of the pdf of a discrete-time Markov chain, $ \psi_t $, given the initial condition $ \psi_0 $.

Start with a simple symmetric tridiagonal matrix

1. Pick some large T and use the initial condition $ \psi_0 = \begin{bmatrix} 1 & 0 & \ldots & 0\end{bmatrix} $
2. Write code to calculate $ \psi_t $ to some $ T $ by iterating the map for each $ t $, i.e.,

$$ \psi_{t+1} = A' \psi_t $$

1. What is the computational order of calculating $ \psi_T $ using this iteration approach $ T < N $?
2. What is the computational order of $ (A')^T = (A' \ldots A') $ and then $ \psi_T = (A')^T \psi_0 $ for $ T < N $?
3. Benchmark calculating $ \psi_T $ with the iterative calculation above as well as the direct $ \psi_T = (A')^T \psi_0 $ to see which is faster. You can take the matrix power with just A_adjoint^T, which uses specialized algorithms faster and more accurately than repeated matrix multiplication (but with the same computational order).
4. Check the same if $ T = 2 N $

Note: The algorithm used in Julia to take matrix powers depends on the matrix structure, as always. In the symmetric case, it can use an eigendecomposition, whereas with a general dense matrix it uses squaring and scaling.

Exercise 2b¶

With the same setup as in Exercise 2a, do an eigendecomposition of A_transpose. That is, use eigen to factor the adjoint $ A' = Q \Lambda Q^{-1} $, where $ Q $ is the matrix of eigenvectors and $ \Lambda $ is the diagonal matrix of eigenvalues. Calculate $ Q^{-1} $ from the results.

Use the factored matrix to calculate the sequence of $ \psi_t = (A')^t \psi_0 $ using the relationship

$$ \psi_t = Q \Lambda^t Q^{-1} \psi_0 $$

where matrix powers of diagonal matrices are simply the element-wise power of each element.

Benchmark the speed of calculating the sequence of $ \psi_t $ up to T = 2N using this method. In principle, the factorization and easy calculation of the power should give you benefits, compared to simply iterating the map as we did in Exercise 2a. Explain why it does or does not, using computational order of each approach.