BAnach Sequential Subspace Optimizer (BASSO)

In the following we assume that you, the reader, has a general familiarity with optimization methods. You should be very well familiar with linear systems of equations and be aware of issues when they are over- oder under-determined, and the basics of inverse problems. You should know what a Tikhonov-Philips functional is, how this can be solved numerically using the Landweber method. You should be familiar with Hilbert spaces and know the principal differences to more general Banach spaces.

If you are not familiar with the above terminology, then we recommend an introductory book to regularization methods such as [Schuster2012]. However, a small introduction is also given in [background].

If you are generally interested in the mathematical background of Banach spaces, Bregman distances and so on, we recommend [Cioranescu1990].

BASSO is a library written in C++ to solve inverse problems in Banach spaces. It is both an interface in C++ and in Python3. Moreover, it comes with a bunch of tools that allow to address certain kind of inverse problems, e.g. inverse problems of ill-conditioned matrices, matrix factorization and certain specific classical inverse problems such as Computerized Tomography. These also may serve as examples how to use the library.

We stress that the library’s interface is written from the viewpoint of the mathematical objects encountered in this problem setup, i.e. there are classes that model vector spaces, normed vector spaces and elements thereof. There are mappings between spaces and functions working on elements from specific spaces. In order to put these to their best use it is recommended to have a background of the mathematical concepts involved.

Therefore, it is not aimed at users who simply wish to call a function to solve a certain type of inverse problem although it comes with a few helper programs that might do the job. At its core it is a library that allows to solve an inverse problem if ones possesses some knowledge of the mathematical underpinning and wishes to see this at work. Extensive care has been taken such that the library automatically prevents mishaps in the tool development such as using vectors from the "wrong" space, use mappings in the wrong way, and so on.

In the following we explain the installation procedure to get BASSO up and running.

As long as no other license statement is given, ThermodynamicAnalyticsToolkit is free for use under the GNU Public License (GPL) Version 2 (see https://www.gnu.org/licenses/gpl-2.0.en.html for full text).

If you encounter any bugs, errors, or would like to submit feature request, please open an issue at GitHub or write to frederik.heber@gmail.com. The authors are especially thankful for any description of all related events prior to occurrence of the error and auxiliary files. More explicitly, the following information is crucial in enabling assistance:

Please mind sensible space restrictions of email attachments.

BASSO mimics all of these objects we have introduced before in its class structure:

These can be considered the basic elements.

Moreover, it contains some more complex classes that comprise the linear system of equations to be solved or the overall problem structure, InverseProblem.

In order to solve this inverse problem of Ax=y we need solvers. Their classes are called LandweberMinimizer, SequentialSubspaceMinimizer or SequentialSubspaceMinimizerNoise (with additional noise in the right-hand side taken into account).

Specific instances of these classes, e.g., the NormedSpace class, are created using the so-called Factory pattern by classes that typically share the same name just with an added Factory, e.g., NormedSpaceFactory. These contain a function create() that is called with certain parameters and returns a pointer to the created instance.

All things associated to an object, e.g., the norm of a space, can be queried from an instance using getters, e.g., NomedSpace.getNorm() for the norm, NormedSpace.getDimension() for the number of degrees of freedom of this space, or NormedSpace.getDualSpace() for its dual space. See [reference.NormedSpace] for details.

The following should give you a rough idea of how to solve an inverse problem

The result will be a structure wherein the approximate solution and a few other instances such as the residual or the number of iterations and the status are contained.

The first four steps can be combined using InverseProblemFactory that takes certain string arguments describing the spaces and the right-hand side and the matrix in Eigen-format.

This is the approach used in the Basso tool that ships with this library, see [examples.basso]

For more complex problems a split feasibility ansatz is used where each constraint is fulfilled in turn in an iterative fashion. In other words, there are two iterations: An inner loop where the inverse problem is solved and an outer loop where the different parts of the split feasibility problem — the inverse problem is one of them — are minimized in turn.

To this end, there is the abstract FeasibilityProblem class, the AuxiliaryConstraintsProblem which combines it with an arbitrary combination of AuxiliaryConstraints such as unity, non-negativity or a logical combination of these. Finally, there is the high-level class SplitFeasibilitySolver that solves inverse problems with these additional auxiliary constraints. These constraints are simply `register()`ed with this class. See [reference.solvers.SplitFeasibilitySolver] for details.

Let us now look at an example where we solve the random matrix toy problem. We create a random matrix and a random right-hand side. We make sure that the right-hand side is actually in the range of the matrix. We construct all instances necessary for the solver and finally solve the inverse problem.

We look at two versions: First, we inspect the "long" version where most of the ingredients are set up by hand. There you will see all of the basic elements that we talked about before. In the second "short" version we will use a lot of convenience objects.

As you probably notice, BASSO is at the moment focused on providing small executables that solve specific (inverse) problems: We have set the options in quite a crude fashion. It is much more elegant to query them from the user and parse them from (argc, argv). Above we have tried to rely on them as little as possible.

Let us look again at the above example. This time we create a CommandLineOptions at the start and use all the convenience objects there are.

You see that all of the basic elements such as setting up of vector spaces and mappings have disappeared. We can still access them through the getters of any of their vectors, see the InverseProblem instance.

These two listings give the basic structure of all example programs in this library. In some cases we combine two minimizations, e.g., we might first project onto the range of the matrix/operator, and then solve the inverse problem without regularization.

In the next section we look at some example programs that have been implemented using the library. We will only discuss the general problem and how the executable is called. Much of it is actually along the same lines as before.

Let us create a vector space and a vector x that is zero initially.

The vector will have 10 components and we are using the $\ell_2$-norm with a power type of 2.

The norm object can be accessed through x.norm, its space by x.space and the dimension can then be obtained by x.space.dim. As mentioned, getters simply are properties. The norm object can be used directly to obtain the norm of x.

Predicates such as SpaceElement::isZero() are available as member functions of the class, i.e., x.isZero().

For creating a LinearMapping, we need two spaces and a matrix.

The matrix needs to be converted from a numpy array to a Eigen::MatrixXd using minieigen. Finally, we create A, the linear mapping.

As a test, let us map an element of space X to space Y.

We see that x has the dimension of X, while the mapped element Ax has the dimension of Y (and is also associated with it). Naturally, the mapping is boring as the matrix is identical to zero in every component. However, we could have used any other numpy array.

The Options class contains all parameters that steer how the iterative solvers work. It is the general parameter storage.

In this case, we instantiate an Options class and set its property algorithm_name to Landweber. This is the name of the iterative solver, here we want to use the Landweber method for solving inverse problems.

For a list of all properties we refer to the respective docstring using help(pyBasso.Options).

Let us now set up an inverse problem. For this, we create a random numpy array that will become our (random) linear operator A. Moreover, we create a random right-hand side y. As before, we need to instantiate two normed vector spaces in order to create a linear mapping using the random matrix.

Again, matrices and vectors are converted from numpy arrays to the Eigen format using minieigen routines.

Now, the random vector is not generally in the range of our random matrix. Therefore, we use a little procedure to have it in its range.

In our case J_p does nothing as we use $\ell_2$-spaces, however in case p is not 2, then we need it.

Finally, we instantiate the inverse problem, set up the options to prepare the solvers and then we solve it.

This will use orthoSESOP to solve for the minimum-norm solution in these Hilbert spaces for up to 1000 steps or if the residuum has a relative change less than $10^{-4}$.

This concludes the discussion of the Python interface of BASSO. All remaining section will now focus again on the C++ parts. However, it should have become clear that this interface lends itself well to solve inverse problems in a straight-forward fashion taking full advantage of BASSO.

The main program which solves any kind of inverse problem given in the form of two files containing the matrix and the right-hand side vector has the same name as the library itself, namely BASSO.

In the following we use it to solve a simple problem which has been investigated in [Schoepfer2006]: Determine the point on a line (in two dimensions) which has minimal distance to the origin in the \ell_p-norm.

We encode the matrix A and the vector y in the files matrix.m and vector.m, respectively.

In this call of Basso we use the Landweber with $\mathrm{X}$ and $\mathrm{Y}$ being \ell_2-spaces. We terminate when the residual is less than 1e-4 and use at most 50 iterations. All information on the iterations is written to an SQLite database file Landweber1.db. This database contains tables data_per_iteration and data_overall tables along with specific views per_iteration and overall, see the class DatabaseManager. The solution, the point closest to the origin on the line described by $A$, is written to straightline1.m.

In the very nice book [Hansen2010] an example of finding the gravitational constant is given. We have implemented this example in the context of BASSO.

This is the problem we have to solve: There is a one-dimensional manifold with "mass" at some distance d to our measurement device. At the surface — the one dimensional manifold where our device sits — we feel the combined effect of the mass on the whole manifold but always diminished by the distance to the respective point mass at f(t). We have information only of the measured gravity field g(s) and now want to deduce the mass distribution f(t). This is our inverse problem: We have the cause. Now we want to know what’s the source of it.

What we need is a equation that relates the two: This is known by the name of Fredholm integral equation of the first kind, see [Hansen2010].

$\int_0^1 K(s,t) f(t) dt = g(s), \quad 0 \leq s \leq 1$,

$K(s,t)$ here is called the "kernel". If the kernel is linear and we discretize the x axis by using a set of equidistant points along it, we obtain a matrix and the equation becomes the typical inverse problems that we know already: $Ax=y$.

$K(s,t) = \frac{d} {\bigl (d^2+(s-t)^2 \bigr)^{\tfrac 3 2}}$.

This equation for the kernel results from knowing that the magnitude of the gravity field along s behaves as $f(t)dt/r^2$, where $r =\sqrt{d^2 + (s-t)^2}$ is the distance between the source point s and the field point t.

This has been encoded in the example Gravity which is called as follows:

The resulting density is stored as a vector in density.m.

In computerized tomography the inverse problem is the reconstruction of the inside of an object from projections. Measurements are for example obtained by passing radiation through a body, whereby their intensity is diminished, proportional to the passed length and density $f(x)$ of the body: [0,1] $^2 \rightarrow R^{+}_{0}$. This decrease is measured over a $a$ angles and $s$ shifts of radiation source and detectors. This measurement matrix is usually called sinogram. Again, we follow [Hansen2010, §7.7].

In (2D) tomography much depends on how the measurement is obtained. Rays might pass through the body in parallel or in fan-like formation. This depends on the setup of the detectors. There might be just a single, point-like detector that is rotated in conjunction with the radiation source around the object on a sphere. The detector might also be spread out, detecting along several small line segments.

In our case we assume single rays, these are suitable parametrized. When the radiation passes through an object, the effect on its intensity is described by the Lambert-Beer law, $dI = I_0 \exp{ f(x) dx}$. Given some initial intensity $I_0$ it is diminished when travelling $dx$ through the body $f(x)$ by $dI$.

This can be written as the famous Radon transform,

$b_i = \int_{-\infty}^{\infty} f(t^i(\tau)) d\tau$,

where i enumerates all rays passed through the body and $t^i$ is the respective parametrization in 2D, $t^i = t^{i,0} + \tau d^i$ with the direction of the ray d and $i=s \cdot a$ rays in total.

The problem is typically discretized in a pixel basis and for setting up the Radon matrix we need to count the length of the ray passing through each of the pixels.

Eventually, we obtain a matrix equation: $b_i = \sum_j a_{ij} x_j$.

Again, this is the inverse problem that we need to solve.

The example program contained with BASSO is called ComputerTomography.

The reconstructed image is written to solution.mat. Changing the norm px changes the reconstruction. It is advised to play around with different values in $(1, \infty)$. Note especially how a value close to 1 for p produces a higher contrast but also more noisy artefacts.

For a problem with a noisy right-hand side we have the following example.

Here, we have the famous Shepp-Logan phantom encoded in phantom.mat which is turned into a right-hand side through $Ax=y$ and perturbed by additional noise of 1% magnitude. The resulting solution is again written to solution.mat. Again, we recommend to play around with px and noise-level.

One last and more complicated program that has been investigated in the course of the BMBF funded project HYPERMATH is the matrix factorization as two alternating inverse problems.

In HYPERMATH we looked at hyperspectral images. These are images that contain not only three colors but thousands of different channels. These images were taken by many different measurement methods: Near-Infrared Spectroscopy (NIR), Matrix Assisted Laser Deposition and Imaging (MALDI), Raman spectroscopy or Spark Optical Emission Spectroscopy (OES). The channels were either truly electromagnetic channels (NIR, Raman), particles masses (MALDI) or elemental abundances (OES). The images per channel consisted of a 2D measurement grid of different objects.

Hyperspectral images are typically very large. $10^6$ to $10^8$ number of pixels and $10^3$ to $10^5$ number of channels are typically encountered, causing up to $10^{12}$ degrees of freedom to be stored. Hence, these are big data.

One crucial problem is therefore storage. The central idea of the HYPERMATH project is to use matrix factorization. The measurements are seen as a large matrix where one dimension is the number of pixels and the other is the number of channels. The underlying assumption is that the measurement is caused by the linear superposition of a few endmembers. For example, in Raman spectroscopy biological tissue is investigated. Spectrographical patterns are caused by the molecules that make up the tissue. There, the endmembers consist of the spectra of each of those individual molecules.

(Non-negative) Matrix factorization will then decouple the superimposed effects of all these endmembers. It splits the measurement matrix A into two factors K and X, where K contains relations from channels to endmembers and X from endmembers to pixels, $A = K \cdot X$.

If the number of reconstructed endmembers is much smaller than the number of pixels or channels, then the storage requirement of both K and X is much smaller than that of the original matrix A.

This problem can be written as two different inverse problems: First, consider each row in X and each row in A. We obtain $A_i = K X_i$, where K is the matrix to formally invert. The other way round, looking at rows, we obtain $A_j = K_j X$. Now, $X^t$ is the matrix to invert.

Next, we perform an alternating least squares minimization, where each inverse problem is solved in turn for a few iterations. This will converge to the two factors we are looking for.

Moreover, this problem can be seen as a split feasibility problem, where each constraint is addressed one by the other. This allows to add further constrains on the factors K and X. For example, row sums might need to be unity or similar.

In the example program MatrixFactorizer this method is implemented.

We only briefly go over the new options:

Eventually, the two low-rank factors are stored in the files spectralmatrix_rank2_nonnegative.m and pixelmatrix_rank2_nonnegative.m.

BASSO comes with several small helper programs that aid in setting up problems, perturbing right hand sides or projecting onto the range of the given matrix.

A finite-dimensional (Banach) space is instantiated using the NormedSpace class by giving its constructor the dimension of the space and a norm object, NormedSpace X(200, norm);`

However, a norm cannot be constructed without a space. Both are tightly interlinked in BASSO (through boost::shared_ptr). Hence, the typical way of setting up a space is by using the NormedSpaceFactory:

Here, we need to specify not only the dimension but also the properties of the norm. We set the type of the norm to be "lp" (or "regularized_l1"), see the static function NormFactory::getMap() for all options.

Moreover, we need to give the space’s parameters. An lp-space needs the value of p and also the power type to create the associated duality mapping. As the number of parameters differ over the different available norms, they are given as a vector of boost::any arguments. Here, we set p to 1.1 and the power type to 2.

Finally, the NormedSpaceFactory::create() function returns the constructed instance of the NormedSpace class. Note that it is wrapped in a boost:shared_ptr and is deconstructed automatically when it is no longer used. In other words, don’t worry about it, just use it like a static instance except you need to dereference it using the * operator on use. Hence, the type has been kept the .._ptr_t ending to remind you that you need to access its members like X->getNorm(), i.e. like a normal pointer.

The SpaceElement class wraps a Eigen::VectorXd and associates it directly with a certain space, see getSpace(). All typical linear algebra routines such as vector product, scalar product, addition, … can be found in this class.

Let us give a few example calls that should be self-explanatory.

The association with the space prevents use in other spaces even if they share the number of degrees of freedom.

A norm basically has the single function to return the "norm" of a given element. The Norm associated with a space makes it a NormedSpace, see getSpace().

Note that certain norms cause spaces to be smooth or not. Hence, the norm can be queried for isSmooth().

Norm only defines the general interface. Specializations are LpNorm, L1Norm, LInfinityNorm that implement the specific norm function to use.

Norms throw an NormIllegalValue_exception when they are called with a parameter outside expected intervals, e.g. an illegal p value such as -1.

Similar to Norm, Mapping defines the general interface of a function that maps an element from one NormedSpace onto another NormedSpace instance, see getSourceSpace() and getTargetSpace(). Therefore, its constructor requires a source and a target space instance.

An example of such a mapping is the LinearMapping. This is a specialization that implements a multiplication with a matrix associated with the mapping.

Its operator() then performs the actual mapping, taking a SpaceElement from the source space X, performing alterations on the element’s components, and returning a SpaceElement in the target space Y.

A more complicated example for using the factory is as follows where we have two matrices whose product forms the actual linear mapping.

Moreover, DualityMapping is another example of a Mapping. There, the constructor only needs a single space as each space knows about its dual space and this mapping always maps into the dual space. In the folder Mappings/Specifics specific implementations of duality mapping, such as the one in lp-spaces, are contained.

Both these classes have specific factory classes, namely LinearMappingFactory and DualityMappingFactory to create the respective instance types. These factories are also used internally, e.g., when constructing the space with its dual space and the associated duality mappings that all have to be interlinked. To this end, these accept "weak" NormedSpace ptr instances.

Finally, each mapping is inherently associated with an adjoint mapping. For LinearMapping this would be the transposed of a (real-value) matrix.

Functions are similar to Mappings and to Norms. They take one or two arguments being SpaceElement instances and produce a real value.

Let us give an example using the BregmanDistance.

Some functions are used solely internally, namely for minimization.

MatrixIO is a namespace wherein several helper functions are gathered that read or write matrix from or to files. Moreover, there are helper functions to print matrices to output streams.

This will parse a vector from the file start_solution.m into the Eigen::VectorXd instance named solution. Here, rhs is used in error messages when parsing, e.., when the file does not exist. Similary, we write the vector into solution.m from the same instance.

When writing SpaceElement`s (instead of `Eigen::VectorXd) we recommend using the convenience class SpaceElementWriter.

The InverseProblem class takes a mapping, two space instances, and a right-hand side. It is simply a convenience struct that keeps all the necessary references at hand, namely it constructs internal shorthand references to all items required for the minimizations methods such as the dual spaces, the adjoint mapping. Moreover, it contains the solution element to be found.

Minimizers form the core of BASSO. They solve the minimizations problems that produce solutions to inverse problems or split feasibility problems.

The general interface is defined in GeneralMinimizer. Its constructor takes all controlling parameters in the form of CommandLineOptions and moreover an InverseProblem instance and a Database whereto iteration-related information is stored. Its operator() expects the inverse problem (as you may change the matrix A and the right), a starting value together with its dual starting value and a possibly known solution to calculate the Bregman distance to it.

The following minimizers are implemented:

These Minimizers may be cleanly instantiated using the MinimizerFactory.

An important point is how the step width is determined. This is especially important for the LandweberMinimizer. There are several different variants implemented. All of which can be easily accessed through the DetermineStepWidthFactory.

The SESOP Minimizers determine the step width automatically through a line search. However, this line search may be tuned to be only inexact by setting appropriate Wolfe conditions.

All Minimizers write information about the iteration to tables of a database. The DatabaseManager controls what values are written and how the tables are set up.

Many parameters that control how for example the minimization is performed are wrapped in a class CommandLineOptions.

Let us briefly discuss what the class looks like and what it can do.

First of all, CommandLineOptions is basically a struct containing all sorts of parameters that control what Banach spaces (and what norms) are employed, what kind of Minimizer, what stopping criteria and so on.

Furthermore, the class has helper functions to quickly parse the state from given command-line options. This explains the nature of its name.

Moreover, the class' state can be stored to a human-readable file using boost::serialization. It is written as in the following example code.

Finally, CommandLineOptions has functions to check the validity of each parameter.

A StoppingCriterion takes several values related to the iteration, namely the time passed, the number of outer iterations, the current residuum and the norm of the right-hand side, and decides whether it is time to stop or not.

There are several primities such as CheckRelativeResiduum or CheckWalltime that may be arbitrarily combined using boolean logic.

It is simplest to employ the factory for producing such a combined stopping criterion. It creates the instance from a simple string as follows.

This will stop either after 100 iterations or when the residuum diverges more than three times (i.e. increases instead of decreases).

See the StoppingCriteriaFactory::StoppingCriteriaFactory() for a list of all known criterion and boolean operators (they are the familar ones from cpp).

The general interface to all solvers is defined in FeasibilityProblem. It is basically simply a functor, i.e. having a single operator() function that executes the solving.

There are some convenience functions such as getName() which returns a human readable name for output purposes. clear() and finish() are used when doing split feasibility problems with multiple iterations.

Then we have several specific solvers which we discuss each in turn: