Matrix mangling for fun and… actually because its not fun and someone else’s solvers are better than yours; they probably have a PhD in it, right?

Solving large systems of equations is not conceptually distinct from optimisation. But, er, there are differences in emphasis between

- finding an optimal vector solution to an equation whose dimension is presumably smaller than your data, and
- finite element methods and PDE wrangling, where the solution might be much larger than your data.

Here I concentrate on the latter.

Clearly artificial neural network require solving what sounds like kind of problem, but they are very very nonlinear, hence kind of a different specialisation. There is some overlap.

Examples applications for solvers can include (especially elliptic) PDEs, inverse/tomography-style image reconstruction, classical optimisation, boring old quadrature, linear transform methods, integral equations, and general sparse-ish problems, various general matrix factorisations…

In many of these domains, there might be structure to exploit, and solvers which can exploit this in some smart way. For example, the description of a discretised PDE usually reduces to a banded matrix and hence a sparse problem.

There are many iterative methods for this kind of system.

Geometric multigrid, algebraic multigrid. Stochastic methods….

fipy is original school:

FiPy is an object oriented, partial differential equation (PDE) solver, written in Python, based on a standard finite volume (FV) approach. … The FiPy framework includes terms for transient diffusion, convection and standard sources, enabling the solution of arbitrary combinations of coupled elliptic, hyperbolic and parabolic PDEs. Currently implemented models include phase field treatments of polycrystalline, dendritic, and electrochemical phase transformations as well as a level set treatment of the electrodeposition process.

pyamg is an algebraic multigrid solver, which is to say, multi-resolution method:

AMG is a multilevel technique for solving large-scale linear systems with optimal or near-optimal efficiency. Unlike geometric multigrid, AMG requires little or no geometric information about the underlying problem and develops a sequence of coarser grids directly from the input matrix. This feature is especially important for problems discretised on unstructured meshes and irregular grids:

A = poisson((500,500), format='csr') # 2D Poisson problem on 500x500 grid ml = ruge_stuben_solver(A) # construct the multigrid hierarchy print ml # print hierarchy information b = rand(A.shape[0]) # pick a random right hand side x = ml.solve(b, tol=1e-10) # solve Ax=b to a tolerance of 1e-8 print "residual norm is", norm(b - A*x) # compute norm of residual vector

See Yousef Saad’s textbooks on algebraic multigrid methods. (free online)

## Refs

- BCFH01
- Brezina, M., Cleary, A., Falgout, R., Henson, V., Jones, J., Manteuffel, T., … Ruge, J. (2001) Algebraic Multigrid Based on Element Interpolation (AMGe).
*SIAM Journal on Scientific Computing*, 22(5), 1570–1592. DOI. - BrHM00
- Briggs, W. L., Henson, V. E., & McCormick, S. F.(2000) A Multigrid Tutorial: Second Edition. . SIAM
- GuWW09
- Guyer, J. E., Wheeler, D., & Warren, J. A.(2009) FiPy: Partial Differential Equations with Python.
*Computing in Science & Engineering*, 11(3), 6–15. DOI. - Stüb01
- Stüben, K. (2001) A review of algebraic multigrid.
*Journal of Computational and Applied Mathematics*, 128(1–2), 281–309. DOI. - TrOS00
- Trottenberg, U., Oosterlee, C. W., & Schuller, A. (2000) Multigrid. . Academic Press
- VaMB96
- Vanek, P., Mandel, J., & Brezina, M. (1996) Algebraic Multigrid By Smoothed Aggregation For Second And Fourth Order Elliptic Problems.
*COMPUTING*, 56, 179–196.