8000 add solve method, tests and update travis instructions. · PythonOptimizers/LLDL.py@d829752 · GitHub
[go: up one dir, main page]
Skip to content

Commit d829752

Browse files
syarrasyarra
committed
add solve method, tests and update travis instructions.
1 parent 28b0e72 commit d829752

File tree

7 files changed

+151
-44
lines changed

7 files changed

+151
-44
lines changed

‎config/.travis.yml

Lines changed: 1 addition & 7 deletions
Original file line numberDiff line numberDiff line change
@@ -13,16 +13,10 @@ before_install:
1313
- virtualenv ~/.venv
1414
- source ~/.venv/bin/activate
1515
- pip install -q numpy
16-
- pip install -q pytest-cov
1716

1817
install:
1918
- cp site.template.cfg site.cfg
2019
- python setup.py install
2120

2221
script:
23-
# - py.test -v --cov-report= --cov=nlp tests/ # when repository goes public
24-
- py.test -v --cov=lldl tests/
25-
26-
# Uncomment when repository goes public
27-
# after_success:
28-
# coveralls
22+
- py.test -v tests/

‎examples/simple_example.py

Lines changed: 6 additions & 18 deletions
Original file line numberDiff line numberDiff line change
@@ -14,7 +14,6 @@
1414

1515
print 'A:'
1616
print A
17-
print A.toarray()
1817

1918
adiag = A.diagonal()
2019
colptrT = np.asarray(T.indptr, dtype=np.int64)
@@ -37,25 +36,14 @@
3736
print u"==== Error : A - LDLᵀ ===="
3837
print ((L * diags(lldl.d, 0) * L.T) - A).toarray()
3938

40-
rhs = A*np.ones([n,1])
39+
print "==== Solving A x = b ===="
40+
rhs = A * np.ones(n)
41+
print "b:"
4142
print rhs
4243

43-
print lldl.colptr
44-
print lldl.rowind
45-
print lldl.lvals
46-
47-
x = lldl.solve(rhs.flatten())
44+
x = lldl.solve(rhs)
4845

46+
print "x:"
4947
print x
50-
#
51-
#x = lldl.solve()
52-
53-
54-
55-
56-
57-
58-
59-
60-
6148

49+
assert np.allclose(x, np.ones(n))

‎examples/simple_example_bis.py

Lines changed: 8 additions & 5 deletions
Original file line numberDiff line numberDiff line change
@@ -13,8 +13,9 @@
1313
print "A:"
1414
print A
1515

16-
solver = CySparseLLDLSolver(A, 48)
17-
(Ls, d, a) = solver.factorize()
16+
n = A.nrow
17+
solver = CySparseLLDLSolver(A, n)
18+
(Ls, d) = solver.factorize()
1819

1920

2021
print 'L:'
@@ -32,11 +33,13 @@
3233
print u"==== Error : A - LDLᵀ ===="
3334
print ((L * D * Lt) - A).to_ll().to_ndarray()
3435

35-
36-
rhs = A*np.ones(48)
36+
print "==== Solving A x = b ===="
37+
rhs = A * np.ones(n)
38+
print "b:"
3739
print rhs
3840

3941
x = solver.solve(rhs)
40-
42+
print "x:"
4143
print x
4244

45+
assert np.allclose(x, np.ones(n))

‎lldl/solver.py

Lines changed: 51 additions & 8 deletions
Original file line numberDiff line numberDiff line change
@@ -12,23 +12,61 @@ class BaseLLDLSolver(object):
1212
"""."""
1313

1414
def __init__(self, n, adiag, colptrT, rowindT, valuesT, memory=5):
15+
"""Instantiate a :class:`BaseLLDLSolver` object.
16+
17+
Perform a limited-memory LDLᵀ factorization of a matrix A
18+
19+
A = L D Lᵀ
20+
where L is a lower triangular matrix and D is a diagonal matrix.
21+
22+
`A` is decomposed as A = T + B + Tᵀ, where T is a strict lower
23+
triangular matrix supplied in CSC (compressed sparse column) format and
24+
B is a diagonal matrix supplied as a one dimensional Numpy array.
25+
26+
The nₖ + p largest elements of L are retained in column k, where nₖ is
27+
the number of nonzero elements in the strict lower triangle of
28+
the k-th column of A and where p ∈ N is a limited-memory factor
29+
specified by the user.
30+
31+
:warning: matrix `A` must be squared.
32+
33+
:parameters:
34+
:n: number of rows (columns) of A
35+
:adiag: diagonal elements of A (as a Numpy array)
36+
:colptrT: pointer to column starts in `rowindT` and `valuesT`
37+
(as a Nympy array)
38+
:rowindT: array of row indices of T (as a Nympy array)
39+
:valuesT: array of non zero elements of T (as a Nympy array)
40+
:memory: limited-memory factor.
41+
"""
1542
self.n = n
1643
self.adiag = adiag.copy()
1744
self.colptrT = colptrT.copy()
1845
self.rowindT = rowindT.copy()
1946
self.valuesT = valuesT.copy()
2047
self.memory = memory
48+
self.factored = False
2149
return
2250

2351
def factorize(self):
24-
(self.colptr, self.rowind, self.lvals, self.d, self.alpha) = sl(self.n, self.n, self.adiag,
25-
self.colptrT, self.rowindT,
26-
self.valuesT, memory=self.memory)
52+
u"""Factorize matrix A as limited-memory LDLᵀ."""
53+
(self.colptr, self.rowind, self.lvals, self.d, self.alpha) = \
54+
sl(self.n, self.n, self.adiag,
55+
self.colptrT, self.rowindT,
56+
self.valuesT, memory=self.memory)
57+
self.factored = True
2758
return
2859

2960
def solve(self, rhs):
61+
u"""Solve A x = b, using the limited-memory LDLᵀ factorization.
62+
63+
:parameters:
64+
:rhs: right-hand side (as a Numpy array).
65+
"""
66+
if not self.factored:
67+
self.factorize()
68+
3069
b = rhs.copy().flatten()
31-
print self.n
3270
dense_lsolve_INT64_FLOAT64(self.n, self.colptr,
3371
self.rowind, self.lvals, b)
3472
z = np.divide(b, self.d)
@@ -41,10 +79,10 @@ class CySparseLLDLSolver(BaseLLDLSolver):
4179
"""Specialized class for CySparse matrices."""
4280

4381
def __init__(self, A, memory=5):
44-
"""Instantiate a :class: `CySparseLLDLSolver`.
82+
"""Instantiate a :class: `CySparseLLDLSolver` object.
4583
4684
:parameters:
47-
:A: the input matrix
85+
:A: the input matrix in CySparse LLmat format.
4886
"""
4987
if not PyLLSparseMatrix_Check(A):
5088
raise TypeError('Input matrix should be a `LLSparseMatrix`')
@@ -55,7 +93,12 @@ def __init__(self, A, memory=5):
5593
A.nrow, A.diag(), colptrT, rowindT, valuesT, memory=memory)
5694

5795
def factorize(self):
58-
"""."""
96+
u"""Factorize matrix A as limited-memory LDLᵀ.
97+
98+
:returns:
99+
:L: L as a llmat matrix
100+
:d: as a Numpy array
101+
"""
59102
super(CySparseLLDLSolver, self).factorize()
60103
nnz = len(self.lvals)
61104
row = np.empty(nnz, dtype=np.int64)
@@ -75,4 +118,4 @@ def factorize(self):
75118
dtype=FLOAT64_T)
76119

77120
L.put_triplet(row, col, val)
78-
return (L, self.d, self.alpha)
121+
return (L, self.d)

‎lldl/src/lldl.cpx

Lines changed: 0 additions & 6 deletions
Original file line numberDiff line numberDiff line change
@@ -49,9 +49,6 @@ def lldl_@index@_@type@(unsigned int n, unsigned int m,
4949
cdef cnp.ndarray[cnp.@index | lower@_t, ndim = 1] rowind = np.zeros(nnzLmax, dtype=np.@index | lower@)
5050
cdef cnp.ndarray[cnp.@index | lower@_t, ndim = 1] colptr = np.zeros(n + 1, dtype=np.@index | lower@)
5151

52-
print colptrT
53-
print rowindT
54-
5552
# Compute the 2-norm of columns of A.
5653
for col in xrange(n):
5754
for k in xrange(colptrT[col], colptrT[col + 1]):
@@ -78,12 +75,9 @@ def lldl_@index@_@type@(unsigned int n, unsigned int m,
7875

7976
success = False
8077
tired = False
81-
print 'colptrT:', colptrT
8278

8379
colptr[:n + 1] = colptrT
8480

85-
print 'colptr:', colptr
86-
8781
while not (success or tired):
8882

8983
# Store the scaled A into L.

‎tests/__init__.py

Whitespace-only changes.

‎tests/test_lldl.py

Lines changed: 85 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,85 @@
1+
# -*- coding: utf-8 -*-
2+
from lldl.solver import BaseLLDLSolver
3+
import numpy as np
4+
from unittest import TestCase
5+
import pytest
6+
try:
7+
from cysparse.sparse.ll_mat import *
8+
import cysparse.common_types.cysparse_types as types
9+
from lldl.solver import CySparseLLDLSolver
10+
except ImportError:
11+
pass
12+
13+
14+
class TestLLDLNympy(TestCase):
15+
16+
def setUp(self):
17+
18+
# Positive definite matrix
19+
self.n = 3
20+
self.indptr = np.array([0, 2, 5, 7])
21+
self.indices = np.array([0, 1, 0, 1, 2, 1, 2])
22+
self.values = np.array([2., -1, -1, 2, -1, -1, 2])
23+
24+
self.colptrT = np.array([0, 1, 2, 2], dtype=np.int64)
25+
self.rowindT = np.array([1, 2], dtype=np.int64)
26+
self.valuesT = np.array([-1., - 1.], dtype=np.float64)
27+
self.adiag = np.array([2., 2., 2.], dtype=np.float64)
28+
29+
def test_init(self):
30+
lldl = BaseLLDLSolver(self.n, self.adiag, self.colptrT,
31+
self.rowindT, self.valuesT, memory=5)
32+
assert lldl.factored is False
33+
34+
def test_factorize(self):
35+
lldl = BaseLLDLSolver(self.n, self.adiag, self.colptrT,
36+
self.rowindT, self.valuesT, memory=5)
37+
lldl.factorize()
38+
assert lldl.factored is True
39+
assert np.array_equal(lldl.colptr, np.array([0, 1, 2, 2],
40+
dtype=np.int64))
41+
assert np.array_equal(lldl.rowind, np.array([1, 2], dtype=np.int64))
42+
assert np.allclose(lldl.lvals, np.array([-0.5, -0.666666667]))
43+
44+
def test_solve(self):
45+
lldl = BaseLLDLSolver(self.n, self.adiag, self.colptrT,
46+
self.rowindT, self.valuesT, memory=5)
47+
rhs = np.array([1, 0, 1], dtype=np.float64)
48+
x = lldl.solve(rhs)
49+
assert np.allclose(x, np.ones(self.n))
50+
51+
52+
class TestLLDLCySparse(TestCase):
53+
54+
def setUp(self):
55+
pytest.importorskip("cysparse.common_types.cysparse_types")
56+
# Positive definite matrix
57+
self.n = 3
58+
self.A = LLSparseMatrix(size=3, itype=types.INT64_T,
59+
dtype=types.FLOAT64_T)
60+
self.Anumpy = np.array([[2., -1., 0.],
61+
[-1., 2., -1.],
62+
[0., -1., 2.]])
63+
self.A[:, :] = self.Anumpy.copy()
64+
65+
def test_init(self):
66+
lldl = CySparseLLDLSolver(self.A, memory=5)
67+
assert lldl.factored is False
68+
69+
def test_factorize(self):
70+
lldl = CySparseLLDLSolver(self.A, memory=5)
71+
(Ls, d) = lldl.factorize()
72+
assert lldl.factored is True
73+
assert PyLLSparseMatrix_Check(Ls) is True
74+
75+
# strict lower triangular matrix is returned.
76+
L = Ls + IdentityLLSparseMatrix(size=Ls.nrow)
77+
Lt = Ls.T + IdentityLLSparseMatrix(size=Ls.nrow)
78+
D = BandLLSparseMatrix(size=Ls.nrow, diag_coeff=[0], numpy_arrays=[d])
79+
assert np.allclose((L * D * Lt).to_ll().to_ndarray(), self.Anumpy)
80+
81+
def test_solve(self):
82+
lldl = CySparseLLDLSolver(self.A, memory=5)
83+
rhs = np.array([1, 0, 1], dtype=np.float64)
84+
x = lldl.solve(rhs)
85+
assert np.allclose(x, np.ones(self.n))

0 commit comments

Comments
 (0)
0