Skip to content

Commit 53a38d8

Browse files
SimonRohouSimonRohou
authored
Merge pull request #378 from damien-masse/codac2_dev
Bugfix for IntFullPivLU::solve
2 parents ce81946 + fa00923 commit 53a38d8

3 files changed

Lines changed: 209 additions & 66 deletions

File tree

‎src/core/matrices/codac2_IntvFullPivLU.cpp‎

Lines changed: 4 additions & 4 deletions
Original file line numberDiff line numberDiff line change
@@ -410,8 +410,9 @@ IntervalMatrix IntvFullPivLU::cokernel() const {
410410
IntervalMatrix IntvFullPivLU::solve(const IntervalMatrix &rhs) const {
411411
assert_release(rhs.rows()==matrixLU_.rows());
412412
IntervalMatrix prhs = _LU.permutationP()*rhs;
413+
Index dim = std::min(matrixLU_.rows(),matrixLU_.cols());
413414
/* inverse L */
414-
for (Index r = 0; r<matrixLU_.rows()-1; r++) {
415+
for (Index r = 0; r<dim; r++) {
415416
for (Index r1=r+1;r1<matrixLU_.rows();r1++) {
416417
prhs.row(r1) -= matrixLU_(r1,r)*prhs.row(r);
417418
}
@@ -428,7 +429,6 @@ IntervalMatrix IntvFullPivLU::solve(const IntervalMatrix &rhs) const {
428429
}
429430
/* then use the diagonal elements */
430431
IntervalMatrix res = IntervalMatrix::Zero(matrixLU_.cols(),rhs.cols());
431-
Index dim = std::min(matrixLU_.cols(),matrixLU_.rows());
432432
for (Index r = dim-1;r>=0;r--) {
433433
for (Index r1=r+1;r1<dim;r1++) {
434434
prhs.row(r) -= matrixLU_(r,r1)*res.row(r1);
@@ -450,8 +450,9 @@ IntervalMatrix IntvFullPivLU::solve(const IntervalMatrix &rhs) const {
450450
void IntvFullPivLU::solve(const IntervalMatrix &rhs, IntervalMatrix &B) const {
451451
assert_release(rhs.rows()==matrixLU_.rows() && B.rows()==matrixLU_.cols());
452452
IntervalMatrix prhs = _LU.permutationP()*rhs;
453+
Index dim = std::min(matrixLU_.cols(),matrixLU_.rows());
453454
/* inverse L */
454-
for (Index r = 0; r<matrixLU_.rows()-1; r++) {
455+
for (Index r = 0; r<dim; r++) {
455456
for (Index r1=r+1;r1<matrixLU_.rows();r1++) {
456457
prhs.row(r1) -= matrixLU_(r1,r)*prhs.row(r);
457458
}
@@ -467,7 +468,6 @@ void IntvFullPivLU::solve(const IntervalMatrix &rhs, IntervalMatrix &B) const {
467468
IntervalMatrix qB = _LU.permutationQ().inverse()*B;
468469
/* then use the diagonal elements */
469470
// IntervalMatrix res = IntervalMatrix::Zero(matrixLU_.cols(),rhs.cols());
470-
Index dim = std::min(matrixLU_.cols(),matrixLU_.rows());
471471
IntervalMatrix U = matrixLU_.triangularView<Eigen::Upper>();
472472
for (Index r = dim-1;r>=0;r--) {
473473
IntervalRow rw=U.row(r);

‎tests/core/matrices/codac2_tests_IntFullPivLU.cpp‎

Lines changed: 65 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -146,3 +146,68 @@ TEST_CASE("IntvFullPivLU")
146146
}
147147

148148
}
149+
150+
TEST_CASE("IntvFullPivLU solve tall matrix")
151+
{
152+
Matrix A({
153+
{ 1, 0 },
154+
{ 0, 1 },
155+
{ 1, 1 },
156+
{ 2,-1 },
157+
{-1, 2 },
158+
});
159+
160+
Matrix X({
161+
{ 2,-1 },
162+
{-1, 3 },
163+
});
164+
165+
Matrix rhs_mid = A*X;
166+
IntervalMatrix rhs(rhs_mid);
167+
168+
IntvFullPivLU lu(A);
169+
170+
SECTION("solve(rhs) with rows >= cols + 2")
171+
{
172+
IntervalMatrix sol = lu.solve(rhs);
173+
174+
CHECK(!sol.is_empty());
175+
CHECK(sol.contains(X));
176+
CHECK((A.template cast<Interval>() * sol).contains(rhs_mid));
177+
}
178+
179+
SECTION("solve(rhs, B) with rows >= cols + 2")
180+
{
181+
IntervalMatrix B = IntervalMatrix::Constant(2,2,{-10.0,10.0});
182+
183+
lu.solve(rhs, B);
184+
185+
CHECK(!B.is_empty());
186+
CHECK(B.contains(X));
187+
CHECK((A.template cast<Interval>() * B).contains(rhs_mid));
188+
}
189+
190+
SECTION("solve(rhs) detects inconsistent rhs")
191+
{
192+
Matrix rhs_bad_mid = rhs_mid;
193+
rhs_bad_mid(4,0) += 1.0; // breaks consistency on the extra row
194+
195+
IntervalMatrix rhs_bad(rhs_bad_mid);
196+
IntervalMatrix sol = lu.solve(rhs_bad);
197+
198+
CHECK(sol.is_empty());
199+
}
200+
201+
SECTION("solve(rhs, B) detects inconsistent rhs")
202+
{
203+
Matrix rhs_bad_mid = rhs_mid;
204+
rhs_bad_mid(4,0) += 1.0; // breaks consistency on the extra row
205+
206+
IntervalMatrix rhs_bad(rhs_bad_mid);
207+
IntervalMatrix B = IntervalMatrix::Constant(2,2,{-10.0,10.0});
208+
209+
lu.solve(rhs_bad, B);
210+
211+
CHECK(B.is_empty());
212+
}
213+
}

‎tests/core/matrices/codac2_tests_IntFullPivLU.py‎

Lines changed: 140 additions & 62 deletions
Original file line numberDiff line numberDiff line change
@@ -13,72 +13,150 @@
1313

1414
class TestIntvFullPivLU(unittest.TestCase):
1515

16+
def test_IntvFullPivLU_1(self):
17+
M = Matrix([
18+
[ 1, -4, 3, 7 ],
19+
[ 2, 1, -4, 6 ],
20+
[ 5, 2, 1 , 9 ],
21+
[ -1, 0, 3, 2 ]
22+
])
1623

17-
# bindings not implemented
18-
def test_IntvFullPivLU_1(self):
19-
M = Matrix([
20-
[ 1, -4, 3, 7 ],
21-
[ 2, 1, -4, 6 ],
22-
[ 5, 2, 1 , 9 ],
23-
[ -1, 0, 3, 2 ]
24-
])
24+
LUdec = IntvFullPivLU(M)
25+
self.assertTrue(LUdec.is_injective()==BoolInterval.TRUE)
26+
self.assertTrue(LUdec.is_surjective()==BoolInterval.TRUE)
27+
self.assertTrue(LUdec.rank()==Interval(4))
28+
self.assertTrue((LUdec.determinant()+602).mag()<=1e-10)
29+
self.assertTrue((LUdec.reconstructed_matrix()-M).norm().ub()<=1e-10)
30+
I1 = LUdec.solve(IntervalMatrix.eye(4,4));
31+
self.assertTrue((I1*M-Matrix.eye(4,4)).norm().ub()<1e-10)
2532

26-
LUdec = IntvFullPivLU(M)
27-
self.assertTrue(LUdec.is_injective()==BoolInterval.TRUE)
28-
self.assertTrue(LUdec.is_surjective()==BoolInterval.TRUE)
29-
self.assertTrue(LUdec.rank()==Interval(4))
30-
self.assertTrue((LUdec.determinant()+602).mag()<=1e-10)
31-
self.assertTrue((LUdec.reconstructed_matrix()-M).norm().ub()<=1e-10)
32-
I1 = LUdec.solve(IntervalMatrix.eye(4,4));
33-
self.assertTrue((I1*M-Matrix.eye(4,4)).norm().ub()<1e-10)
33+
def test_IntvFullPivLU_2(self):
34+
M = Matrix([
35+
[ 1, -4, 6, 7 ],
36+
[ 2, 1, 3, 6 ],
37+
[ 5, 2, 8 , 9 ],
38+
[ -1, 0, -2, 2 ]
39+
])
3440

35-
def test_IntvFullPivLU_2(self):
36-
M = Matrix([
37-
[ 1, -4, 6, 7 ],
38-
[ 2, 1, 3, 6 ],
39-
[ 5, 2, 8 , 9 ],
40-
[ -1, 0, -2, 2 ]
41-
])
41+
LUdec = IntvFullPivLU(M)
42+
self.assertTrue(LUdec.is_injective()==BoolInterval.UNKNOWN)
43+
self.assertTrue(LUdec.is_surjective()==BoolInterval.UNKNOWN)
44+
self.assertTrue(LUdec.rank()==Interval([3,4]))
45+
self.assertTrue((LUdec.determinant()).mag()<=1e-10)
46+
self.assertTrue((LUdec.reconstructed_matrix()-M).norm().ub()<=1e-10)
47+
K = LUdec.kernel()
48+
self.assertTrue(K.cols()==1)
49+
self.assertTrue((M*K).norm().ub()<1e-10)
50+
coK = LUdec.cokernel()
51+
self.assertTrue(coK.rows()==1)
52+
self.assertTrue((coK*M).norm().ub()<1e-10)
53+
Im = LUdec.image(M)
54+
self.assertTrue(Im.cols()==3)
55+
coIm = LUdec.coimage(M)
56+
self.assertTrue(coIm.rows()==3)
4257

43-
LUdec = IntvFullPivLU(M)
44-
self.assertTrue(LUdec.is_injective()==BoolInterval.UNKNOWN)
45-
self.assertTrue(LUdec.is_surjective()==BoolInterval.UNKNOWN)
46-
self.assertTrue(LUdec.rank()==Interval([3,4]))
47-
self.assertTrue((LUdec.determinant()).mag()<=1e-10)
48-
self.assertTrue((LUdec.reconstructed_matrix()-M).norm().ub()<=1e-10)
49-
K = LUdec.kernel()
50-
self.assertTrue(K.cols()==1)
51-
self.assertTrue((M*K).norm().ub()<1e-10)
52-
coK = LUdec.cokernel()
53-
self.assertTrue(coK.rows()==1)
54-
self.assertTrue((coK*M).norm().ub()<1e-10)
55-
Im = LUdec.image(M)
56-
self.assertTrue(Im.cols()==3)
57-
coIm = LUdec.coimage(M)
58-
self.assertTrue(coIm.rows()==3)
58+
def test_IntvFullPivLU_3(self):
59+
M = Matrix([
60+
[ 1, -4, 6, 7, 6 ],
61+
[ 2, 1, 3, 6, -2 ],
62+
[ 5, 2, 8 , 9, -1 ]
63+
])
64+
LUdec = IntvFullPivLU(M)
65+
self.assertTrue(LUdec.is_injective()==BoolInterval.FALSE)
66+
self.assertTrue(LUdec.is_surjective()==BoolInterval.TRUE)
67+
self.assertTrue(LUdec.rank()==Interval(3))
68+
self.assertTrue((LUdec.reconstructed_matrix()-M).norm().ub()<=1e-10)
69+
K = LUdec.kernel()
70+
self.assertTrue(K.cols()==2)
71+
self.assertTrue((M*K).norm().ub()<1e-10)
72+
Im = LUdec.image(M)
73+
self.assertTrue(Im.cols()==3)
74+
I1 = LUdec.solve(IntervalMatrix.eye(3,3));
75+
self.assertTrue((M*I1-Matrix.eye(3,3)).norm().ub()<1e-10)
76+
A = IntervalMatrix([ [1], [[-20,20]], [2], [[-20,20]], [[-20,20]] ])
77+
B = IntervalMatrix([ [2.0], [1.0], [4.0] ])
78+
LUdec.solve(B,A)
79+
self.assertTrue((M*A-B).norm().ub()<=1e-10)
5980

60-
def test_IntvFullPivLU_3(self):
61-
M = Matrix([
62-
[ 1, -4, 6, 7, 6 ],
63-
[ 2, 1, 3, 6, -2 ],
64-
[ 5, 2, 8 , 9, -1 ]
65-
])
66-
LUdec = IntvFullPivLU(M)
67-
self.assertTrue(LUdec.is_injective()==BoolInterval.FALSE)
68-
self.assertTrue(LUdec.is_surjective()==BoolInterval.TRUE)
69-
self.assertTrue(LUdec.rank()==Interval(3))
70-
self.assertTrue((LUdec.reconstructed_matrix()-M).norm().ub()<=1e-10)
71-
K = LUdec.kernel()
72-
self.assertTrue(K.cols()==2)
73-
self.assertTrue((M*K).norm().ub()<1e-10)
74-
Im = LUdec.image(M)
75-
self.assertTrue(Im.cols()==3)
76-
I1 = LUdec.solve(IntervalMatrix.eye(3,3));
77-
self.assertTrue((M*I1-Matrix.eye(3,3)).norm().ub()<1e-10)
78-
A = IntervalMatrix([ [1], [[-20,20]], [2], [[-20,20]], [[-20,20]] ])
79-
B = IntervalMatrix([ [2.0], [1.0], [4.0] ])
80-
LUdec.solve(B,A)
81-
self.assertTrue((M*A-B).norm().ub()<=1e-10)
81+
def test_IntvFullPivLU_solve_tall_matrix(self):
8282

83+
A = Matrix([
84+
[ 1, 0 ],
85+
[ 0, 1 ],
86+
[ 1, 1 ],
87+
[ 2,-1 ],
88+
[-1, 2 ],
89+
])
90+
91+
X = Matrix([
92+
[ 2,-1 ],
93+
[-1, 3 ],
94+
])
95+
96+
rhs_mid = A*X
97+
rhs = IntervalMatrix(rhs_mid)
98+
99+
lu = IntvFullPivLU(A)
100+
101+
sol = lu.solve(rhs)
102+
103+
self.assertTrue(not sol.is_empty())
104+
self.assertTrue(sol.contains(X))
105+
self.assertTrue((A*sol).contains(rhs_mid))
106+
107+
def test_IntvFullPivLU_solve_tall_matrix_with_box(self):
108+
109+
A = Matrix([
110+
[ 1, 0 ],
111+
[ 0, 1 ],
112+
[ 1, 1 ],
113+
[ 2,-1 ],
114+
[-1, 2 ],
115+
])
116+
117+
X = Matrix([
118+
[ 2,-1 ],
119+
[-1, 3 ],
120+
])
121+
122+
rhs_mid = A*X
123+
rhs = IntervalMatrix(rhs_mid)
124+
125+
lu = IntvFullPivLU(A)
126+
B = IntervalMatrix.constant(2,2,[-10,10])
127+
128+
lu.solve(rhs, B)
129+
130+
self.assertTrue(not B.is_empty())
131+
self.assertTrue(B.contains(X))
132+
self.assertTrue((A*B).contains(rhs_mid))
133+
134+
def test_IntvFullPivLU_solve_tall_matrix_inconsistent_rhs(self):
135+
136+
A = Matrix([
137+
[ 1, 0 ],
138+
[ 0, 1 ],
139+
[ 1, 1 ],
140+
[ 2,-1 ],
141+
[-1, 2 ],
142+
])
143+
144+
rhs_bad = IntervalMatrix([
145+
[[ 2, 2],[-1,-1]],
146+
[[-1,-1],[ 3, 3]],
147+
[[ 1, 1],[ 2, 2]],
148+
[[ 5, 5],[-5,-5]],
149+
[[-3,-3],[ 7, 7]], # should be [-4, 7] for consistency
150+
])
151+
152+
lu = IntvFullPivLU(A)
153+
154+
sol = lu.solve(rhs_bad)
155+
self.assertTrue(sol.is_empty())
156+
157+
B = IntervalMatrix.constant(2,2,[-10,10])
158+
lu.solve(rhs_bad, B)
159+
self.assertTrue(B.is_empty())
160+
83161
if __name__ == '__main__':
84-
unittest.main()
162+
unittest.main()

0 commit comments

Comments
 (0)