Numerical Linear Algebra and Optimization

This is the fourth post in our series on Krylov subspaces. The previous ones (i.e Arnoldi Iterations and Lanczos Algorithm were mostly focused on eigenvalue and eigenvector computations. In this post we will have a look at solving strategies for linear systems of equations. In particular, we will derive the famous conjugate gradients algorithm. The conjugate gradients method is one of the most popular solvers for linear systems of equations. There exist countless variations and extensions.

This is the third post in my series on Krylov subspaces. The first post is here and the second one is here. The Lanczos Algorithm In this post we cover the Lanczos algorithm that gives you eigenvalues and eigenvectors of symmetric matrices. The Lanczos algorithm is named after Cornelius Lanczos, who was quite influential with his research. He also proposed an interpolation method and a method to approximate the gamma function.

This is the second post in my series on Krylov subspaces. The first post is here and the third one is here. Arnoldi Iterations Arnoldi iterations is an algorithm to find eigenvalues and eigenvectors of general matrices. The algorithm constructs orthonormal bases of suitable Krylov subspaces and extracts the eigendata from the basis representation. The method is named after the American engineer Walter Edwin Arnoldi and was published in 1951 [1].

This is the first post in a (planned) series on Krylov subspaces, projection processes, and related algorithms. We already discussed projection processes when talking about the Bregman algorithm and we will see that the Krylov (sub-)spaces will be generated by a set of vectors that are not necessarily orthogonal. Getting an orthogonal basis for these spaces can for example be done by applying Gram-Schmidt. Forthcoming posts will discuss further interesting applications such as the Arnoldi iteration or the Lanczos method.

In the previous post in our series on the Bregman algorithm we discussed how to solve convex optimization problems. In this post we want to give reference to some variations and extensions of the Bregman algorithm. First of all, there’s a website related to the Split Bregman Algorithm, which we analysed here. The website provides further documentation and reference implementations in Matlab. There is also a nice paper [1] that discusses its relationship to other frameworks in convex optimisation.

In a previous post we discussed how to solve constrained optimization problems by using the Bregman algorithm. Here we want to extend the approach unconstrained problems. Let’s start simple. Assume we want to minimize a convex and smooth function $f\colon\mathbb{R}^{n}\to\mathbb{R}$. Since $\mathbb{R}^n$ is a convex set, we can apply the Bregman algorithm from our previous post and get the following iterative scheme. $$\begin{align} x^{(n+1)} &= \arg\min_z\left\lbrace D(z, x^{(n)})\right\rbrace \\ &= \arg\min_z\left\lbrace f(z) - \left\langle \nabla f(x^{(n)}), z \right\rangle \right\rbrace \end{align}$$

In a previous post we discussed how to find a common point in a family of convex sets by using the Bregman algorithm. Actually the algorithm is capable of more. We can use it to solve constrained optimization problems. Since there is no change on the underlying algorithm, we will focus on an example here. Assume we are given a quadratic form $f(x) = x^\top Q x$ with a symmetric positive definite matrix $Q$.

In the 1960s Lev Meerovich Bregman developed an optimization algorithm [1] which became rather popular beginning of 2000s. It’s not my intention to present the proofs for all the algorithmic finesse, but rather the general ideas why it is so appealing. Some ramifications and new algorithms, that were derived from the original Bregman algorithm (e.g. [2, 3, 4]), will be presented in subsequent posts. Finding a common point in a set of convex sets Let’s assume we are given a family of convex sets $A_i$, with $i=1, \ldots, n$.

We discuss the necessary changes to our QR decomposition algorithms to handle matrices which do not have full rank.

We developed a QR decomposition algorithm, based on the orthogonalisation process of Gram-Schmidt in a series of posts here, here, here, and here. Let’s have a look how good this algorithm performs against built-in implementations from julia and other programming languages. Test criteria For a given matrix $A$ we get a decomposition $A=QR$ where $Q$ is orthogonal and where $R$ is upper triangular. If $A$ does not have full column rank, then some columns in our matrix $Q$ will be 0 (at least in our implementation).

Problem Formulation We already discussed QR decompositions and showed that using the modified formulation of Gram-Schmidt significantly improves the accuracy of the results. However, there is still an error of about $10^3 M_\varepsilon$ (where $M_\varepsilon$ is the machine epsilon) when using the modified Gram Schmidt as base algorithm for the orthogonalisation. These errors are due to cancellations due to the limited precision in our floating point representation. Reorthogonalisations We consider the following idea: each time we compute the next orthogonal vector, we check whether cancellations occurred and whether our result might be inaccurate.

We consider the necessary changes to the Gram-Schmidt orthogonalisation to obtain a QR Decomposition

We compare the accuracy of the classical Gram-Schmidt algorithm to the modified Gram-Schmidt algorithm.