hub.darcs.net :: thielema -> lapack -> files

Installation

Before installing the Haskell bindings you need to install the BLAS and LAPACK packages. Please note, that additionally to the reference implementation in FORTRAN 77 there are alternative optimized implementations like OpenBLAS, ATLAS, Intel MKL.

Debian, Ubuntu

sudo apt-get install libblas-dev liblapack-dev

MacOS

You may install pkgconfig and LAPACK via MacPorts:

sudo port install pkgconfig lapack

However, the pkg-config files for LAPACK seem to be installed in a non-standard location. You must make them visible to pkg-config by

export PKG_CONFIG_PATH=/opt/local/lib/lapack/pkgconfig:$PKG_CONFIG_PATH

You may set the search PATH permanently by adding the above command line to your $HOME/.profile file.

Alternatively, a solution for all users of your machine would be to set symbolic links:

sudo ln -s /opt/local/lib/lapack/pkgconfig/blas.pc /opt/local/lib/pkgconfig/blas.pc
sudo ln -s /opt/local/lib/lapack/pkgconfig/lapack.pc /opt/local/lib/pkgconfig/lapack.pc

Introduction

Here is a small example for constructing and formatting matrices:

Prelude> import qualified Numeric.LAPACK.Matrix as Matrix
Prelude Matrix> import Numeric.LAPACK.Format as Fmt ((##))
Prelude Matrix Fmt> let a = Matrix.fromList (Matrix.shapeInt 3) (Matrix.shapeInt 4) [(0::Float)..]
Prelude Matrix Fmt> a ## "%.4f"
 0.0000 1.0000  2.0000  3.0000
 4.0000 5.0000  6.0000  7.0000
 8.0000 9.0000 10.0000 11.0000
Prelude Matrix Fmt> import qualified Numeric.LAPACK.Matrix.Shape as MatrixShape
Prelude Matrix Fmt MatrixShape> import qualified Numeric.LAPACK.Matrix.Triangular as Triangular
Prelude Matrix Fmt MatrixShape Triangular> let u = Triangular.upperFromList MatrixShape.RowMajor (Matrix.shapeInt 4) [(0::Float)..]
Prelude Matrix Fmt MatrixShape Triangular> (u, Triangular.transpose u) ## "%.4f"
 0.0000 1.0000 2.0000 3.0000
        4.0000 5.0000 6.0000
               7.0000 8.0000
                      9.0000

 0.0000
 1.0000 4.0000
 2.0000 5.0000 7.0000
 3.0000 6.0000 8.0000 9.0000

You may find a more complex introductory example at: http://code.henning-thielemann.de/bob2019/main.pdf

Formatting

We do not try to do fancy formatting for the Show instance. The Show instances of matrices emit plain valid Haskell code in one line where all numbers are printed in full precision. If matrices are part of larger Haskell data structures this kind of formatting works best. For human-friendly formatting in GHCi you need to add something like ## "%.4f" after a matrix or vector expression. It means: Print all numbers in fixed point representation using four digits for the fractional part. You can use the formatting placeholders provided by printf. The matrices have Hyper instances for easy usage in HyperHaskell.

Formatting is based on the Output type class that currently supports output as Text boxes for GHCi and HTML for HyperHaskell.

You may tell GHCi to use Format class instead of Show:

Fmt> let lapackPrint x = x##"%.3f"
Fmt> :set -interactive-print lapackPrint

You may permanently configure this one in your local .ghci file. If you want to display values via Show class, you can always fall back by:

Fmt> print "Hello"

Matrix vs. Vector

Vectors are Storable.Arrays from the comfort-array package. An array can have a fancy shape like a shape defined by an enumeration type, the shape of two appended arrays, a shape that is compatible to a Haskell container type, a rectangular or triangular shape.

All operations check dynamically whether corresponding shapes match structurally. E.g. a|||b === c|||d composes a matrix from four quadrants a, b, c, d. It is not enough that a|||b and c|||d have the same width, but the widths of a and c as well as b and d must match. The type variables for shapes show which dimensions must be compatible. We recommend to use type variables for the shapes as long as possible because this increases type safety even if you eventually use only ShapeInt. If you use statically sized shapes you get static size checks.

A matrix can have any of these shapes as height and as width. E.g. it is no problem to define a matrix that maps a triangular shaped array to a rectangular shaped one. There are actually practical applications to such matrices. A matrix can be treated as vector, but there are limitations. E.g. if you scale a Hermitian matrix by a complex factor it will in general be no longer Hermitian. Another problem: Two equally sized rectangular matrices may differ in the element order (row major vs. column major). You cannot simply add them by adding the flattened arrays element-wise. Thus if you want to perform vector operations on a matrix the package requires you to "unpack" a matrix to a vector using Matrix.Array.toVector. This conversion is almost a no-op and preserves most of the shape information. The reverse operation is Matrix.Array.fromVector.

There are more matrix types that are not based on a single array. E.g. we provide a symbolic inverse, a scaling matrix, a permutation matrix. We also support arrays that represent factors of a matrix factorization. You obtain these by LU and QR decompositions. You can extract the matrix factors of it, but you can also multiply the factors to other matrices or use the decompositions for solving simultaneous linear equations.

Matrix type parameters

LAPACK supports a variety of special matrix types, e.g. dense, banded, triangular, symmetric, Hermitian matrices and our Haskell interface supports them, too. There are two layers: The low level layer addresses how matrices are stored for LAPACK. Matrices and vectors are stored in the Array type from comfort-array:Data.Array.Comfort.Storable using shape types from Numeric.LAPACK.Matrix.Layout. The high level layer provides a matrix type in Numeric.LAPACK.Matrix.Array with mathematically relevant type parameters. The matrix type is:

ArrayMatrix pack property lower upper meas vert horiz height width a

The type parameters are from right to left:

• a - the element type

• height and width are the vertical and horizontal shapes of the matrix

• meas vert horiz form a group with following possible assignments:

  meas	vert	horiz	name	condition
  Shape	Small	Small	Square matrix	height == width
  Size	Small	Small	Liberal square	size height == size width
  Size	Big	Small	Tall matrix	size height >= size width
  Size	Small	Big	Wide matrix	size height <= size width
  Size	Big	Big	General matrix	arbitrary height and width

  Think of meas as the measurement that goes into the relation of dimensions. You can either compare shapes (meas ~ Shape) or their sizes (meas ~ Size). For vert and horiz the possible values mean:

  • Small: The corresponding dimension is equal to the minimum of height and width.

  • Big: The corresponding dimension has no further restrictions, but it is of course at least the minimum of height and width.

  The names Small and Big fit best to tall and wide matrices. The remaining combinations Small Small for Square and Big Big for General appear to be arbitrary, but they help to e.g. treat square and tall matrices the same way, where sensible. Turning Shape into Size or Small into Big relaxes a dimension relation.

• lower upper count the numbers of non-zero off-diagonals.

  Of course, stored off-diagonals can consist entirely of zeros. Thus more precisely we have to say, that lower and upper tell that all values outside the numbered bands are zero.

  lower and upper can be:

  • Filled - no restriction on the number of off-diagonals.

  • Bands n, where n is a natural number unarily encoded in types.

  Empty is a synonym for Bands U0.

• property can be

  • Arbitrary - this type does not make any further promises about the matrix elements

  • Unit - matrix is triangular with a unit diagonal

  It can be used for matrices that always have a unit diagonal by construction. This property is preserved by some operations and enables optimizations by LAPACK. Solving with a unit triangular matrix does not require division and thus cannot fail due to division by zero.

  • Symmetric - matrix is symmetric

  • Hermitian - matrix is Hermitian (also supported for real elements)

  The internal Hermitian property also has three type tags neg zero pos to restrict the range of values of bilinear-forms. This way you can denote positive definiteness and semidefiniteness.

• pack can be Packed or Unpacked.

Unpacked means that the full matrix bounded by height and width is stored. Packed format is supported for triangular, symmetric, Hermitian and banded matrices.

For banded matrices you should always prefer the packed format. For triangular, symmetric and Hermitian matrices LAPACK does not always support packed format natively and our Haskell interface converts forth and back silently. I also think that unpacked triangular formats enjoy better support by vectorized block algorithms. Thus, the unpacked triangular format may be better for performance although it requires twice as much space as the packed format. The packed triangular format however is still the default format for conversion to and from lists, because this prevents the user from declaring non-zero values in the empty area of a triangular matrix.

Let us examine some examples:

• Diagonal matrix:

  ArrayMatrix Packed Arbitrary Empty Empty Shape Small Small sh sh a

• Packed upper triangular matrix:

  ArrayMatrix Packed Arbitrary Empty Filled Shape Small Small sh sh a

• Unpacked unit lower triangular matrix:

  ArrayMatrix Unpacked Unit Filled Empty Shape Small Small sh sh a

• Complex-valued symmetric matrix:

  ArrayMatrix Packed Symmetric Filled Filled Shape Small Small sh sh (Complex a)

• Tall banded matrix:

  ArrayMatrix Packed Arbitrary (Bands sub) (Bands super) Size Big Small height width a

The type tags have a mathematical meaning and this pays off for operations on matrices. E.g. matrix multiplication adds off-diagonals. Matrix inversion fills non-zero triangular matrix parts. The five supported relations for matrix dimensions are transitive, and thus matrix multiplication maintains a dimension relation, e.g. tall times tall is tall.

Please note, that not all type parameter combinations are supported. Some restrictions are dictated by mathematics, e.g. Hermitian matrices must always be square, matrices with unit diagonal must always be triangular and so on. Some combinations are simply not supported, because they do not add value. E.g. a (square) diagonal matrix is always symmetric but we allow Symmetric only together with Filled. Forbidden combinations are often not prevented at the type level, but you will not be able to construct a matrix of a forbidden type.

Matrix type synonyms

• General: A dense rectangular matrix with no restrictions to the sizes.

• Tall: A dense rectangular matrix with size height >= size width.

• Wide: A dense rectangular matrix with size height <= size width.

• Square: A dense square matrix with height == width. If you miss a function in Numeric.LAPACK.Matrix it is likely that it is actually a function only defined for Square matrices. In this case, please check Numeric.LAPACK.Matrix.Square.

• LiberalSquare: A dense square matrix with size height == size width.

• Full: A dense matrix that might be either General, Tall, Wide, Square, LiberalSquare.

• Triangular: An upper or lower triangular matrix with height == width.

• Symmetric, Hermitian: If you miss a function in Symmetric it might be defined in Hermitian instead, where Symmetric and Hermitian are the same for real-valued matrices.

• Mosaic: A Quadratic matrix that is either Triangular, Symmetric, Hermitian.

• Quadratic: A square matrix with height == width that might be Full, Triangular, Symmetric, Hermitian, Banded, BandedHermitian.

Infix operators

The package provides fancy infix operators like #*| and \*#. They symbolize both operands and operations. E.g. in #*| the hash means Matrix, the star means Multiplication and the bar means Column Vector.

Possible operations are:

• a_*_b - a multiplied by b

• a_/_b - a multiplied by inverse b

• a_\_b - inverse a multiplied by b

Possible operands are:

• # - a matrix that is generalized through a type class

• ## - a full matrix

• \ - a diagonal matrix represented by a Vector

• - - a row vector

• | - a column vector

• . - a scalar

For multiplication of equally shaped matrices we also provide instances of Semigroup.<>.

Precedence of the operators is chosen analogously to plain * and /. Associativity is chosen such that the same operator can be applied multiple times without parentheses. But sometimes this may mean that you have to mix left and right associative operators, and thus you may still need parentheses.

Type errors

You might encounter cryptic type errors that refer to the encoding of particular matrix types via matrix type parameters.

E.g. the error

Couldn't match type `Numeric.LAPACK.Matrix.Extent.Big`
               with `Numeric.LAPACK.Matrix.Extent.Small`

may mean that you passed Square where General or Tall was expected. You may solve the problem with a function like Square.toFull or Square.fromFull.

The error

Couldn't match type `Type.Data.Bool.False`
               with `Type.Data.Bool.True`

most likely refers to non-matching definiteness warranties in a Hermitian matrix. You may try a function like Hermitian.assureFullRank or Hermitian.relaxIndefinite to fix the issue.

Tips and Tricks

Solving simultaneous linear equations

LAPACK provides many solvers, for each type of matrix at least one, sometimes more. The Haskell lapack package should provide interfaces to all of them. However, there are two small problems: First, all solvers solve for many right-hand sides at once. Many right-hand side vectors must be assembled as columns of a matrix. However, in most applications you have only one right-hand side. You can adapt a matrix-rhs solver to single vectors like so: Matrix.unliftColumn Layout.ColumnMajor (solve matrix) rhs. This is cumbersome. Second, the infix operator #\| does the unlifting for you: matrix #\| rhs. It also chooses the default solver for the kind of matrix. This choice might not be appropriate for you, since e.g. Gauss elimination can be numerically imprecise. You can choose a different solver for #\| by doing the matrix decomposition manually. This looks like so: Householder.fromMatrix matrix #\| rhs.