1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
|
Solves the linear system JtJ x = b using 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.
As many vectors b as we'd like may be given at one time (in rows of bt). The
dimensions of the returned array xt will match the dimensions of the given array
bt.
Broadcasting is supported: any leading dimensions will be processed correctly,
as long as bt has shape (..., Nstate)
This function 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:
- bt must be C-contiguous (the normal numpy order)
- bt must contain 64-bit floating-point values (dtype=float)
This function is now able to pass different values of "sys" to the internal
cholmod_solve2() call. This is specified with the "mode" argument. By default,
we use CHOLMOD_A, which is the default behavior: we solve JtJ x = b. All the
other modes supported by CHOLMOD are supported. From cholmod.h:
CHOLMOD_A: solve Ax=b
CHOLMOD_LDLt: solve LDL'x=b
CHOLMOD_LD: solve LDx=b
CHOLMOD_DLt: solve DL'x=b
CHOLMOD_L: solve Lx=b
CHOLMOD_Lt: solve L'x=b
CHOLMOD_D: solve Dx=b
CHOLMOD_P: permute x=Px
CHOLMOD_Pt: permute x=P'x
See the CHOLMOD documentation and source for details.
ARGUMENTS
- bt: a numpy array of shape (..., Nstate). This array must be C-contiguous and
it must have dtype=float
- sys: optional string, defaulting to "A": solve JtJ x = b. Selects the specific
problem being solved; see the description above. The value passed to "sys"
should be the string with or without the "CHOLMOD_" prefix
RETURNED VALUE
The transpose of the solution array x, in a numpy array of the same shape as the
input bt
|
