CHOLMOD_factorization.docstring

A basic Python interface to CHOLMOD

SYNOPSIS

    from scipy.sparse import csr_matrix

    indptr  = np.array([0, 2, 3, 6, 8])
    indices = np.array([0, 2, 2, 0, 1, 2, 1, 2])
    data    = np.array([1, 2, 3, 4, 5, 6, 7, 8], dtype=float)

    Jsparse = csr_matrix((data, indices, indptr))
    Jdense  = Jsparse.toarray()
    print(Jdense)
    ===> [[1. 0. 2.] 
          [0. 0. 3.] 
          [4. 5. 6.] 
          [0. 7. 8.]]

    bt = np.array(((1., 5., 3.), (2., -2., -8)))
    print(nps.transpose(bt))
    ===> [[ 1.  2.] 
          [ 5. -2.] 
          [ 3. -8.]]

    F  = mrcal.CHOLMOD_factorization(Jsparse)
    xt = F.solve_xt_JtJ_bt(bt)
    print(nps.transpose(xt))
    ===> [[ 0.02199662  0.33953751] 
          [ 0.31725888  0.46982516] 
          [-0.21996616 -0.50648618]]

    print(nps.matmult(nps.transpose(Jdense), Jdense, nps.transpose(xt)))
    ===> [[ 1.  2.] 
          [ 5. -2.] 
          [ 3. -8.]]

The core of the mrcal optimizer is a sparse linear least squares solver using
CHOLMOD to solve a large, sparse linear system. CHOLMOD is a C library, but it
is sometimes useful to invoke it from Python.

The CHOLMOD_factorization class factors a matrix JtJ, and this method uses that
factorization to efficiently solve the linear equation JtJ x = b. The usual
linear algebra conventions refer to column vectors, but numpy generally deals
with row vectors, so I talk about solving the equivalent transposed problem: xt
JtJ = bt. The difference is purely notational.

The class takes a sparse array J as an argument in __init__(). J is optional,
but there's no way in Python to pass it later, so from Python you should always
pass J. This is optional for internal initialization from C code.

J must be given as an instance of scipy.sparse.csr_matrix. csr is a row-major
sparse representation. CHOLMOD wants column-major matrices, so it see this
matrix J as a transpose: the CHOLMOD documentation refers to this as "At". And
the CHOLMOD documentation talks about factoring AAt, while I talk about
factoring JtJ. These are the same thing.

The factorization of JtJ happens in __init__(), and we use this factorization
later (as many times as we want) to solve JtJ x = b by calling
solve_xt_JtJ_bt().

This class carefully checks its input for validity, but makes no effort to be
flexible: anything that doesn't look right will result in an exception.
Specifically:

- J.data, J.indices, J.indptr must all be numpy arrays

- J.data, J.indices, J.indptr must all have exactly one dimension

- J.data, J.indices, J.indptr must all be C-contiguous (the normal numpy order)

- J.data must hold 64-bit floating-point values (dtype=float)

- J.indices, J.indptr must hold 32-bit integers (dtype=np.int32)

ARGUMENTS

The __init__() function takes

- J: a sparse array in a scipy.sparse.csr_matrix object