/*
* Licensed to the Apache Software Foundation (ASF) under one or more
* contributor license agreements. See the NOTICE file distributed with
* this work for additional information regarding copyright ownership.
* The ASF licenses this file to You under the Apache License, Version 2.0
* (the "License"); you may not use this file except in compliance with
* the License. You may obtain a copy of the License at
*
* http://www.apache.org/licenses/LICENSE-2.0
*
* Unless required by applicable law or agreed to in writing, software
* distributed under the License is distributed on an "AS IS" BASIS,
* WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
* See the License for the specific language governing permissions and
* limitations under the License.
*/
package org.apache.commons.math4.linear;
import java.util.Arrays;
import org.apache.commons.math4.exception.DimensionMismatchException;
import org.apache.commons.math4.linear.Array2DRowRealMatrix;
import org.apache.commons.math4.linear.ArrayRealVector;
import org.apache.commons.math4.linear.DecompositionSolver;
import org.apache.commons.math4.linear.IterativeLinearSolver;
import org.apache.commons.math4.linear.IterativeLinearSolverEvent;
import org.apache.commons.math4.linear.JacobiPreconditioner;
import org.apache.commons.math4.linear.LUDecomposition;
import org.apache.commons.math4.linear.NonPositiveDefiniteOperatorException;
import org.apache.commons.math4.linear.NonSelfAdjointOperatorException;
import org.apache.commons.math4.linear.NonSquareOperatorException;
import org.apache.commons.math4.linear.PreconditionedIterativeLinearSolver;
import org.apache.commons.math4.linear.RealLinearOperator;
import org.apache.commons.math4.linear.RealVector;
import org.apache.commons.math4.linear.SymmLQ;
import org.apache.commons.math4.util.FastMath;
import org.apache.commons.math4.util.IterationEvent;
import org.apache.commons.math4.util.IterationListener;
import org.junit.Assert;
import org.junit.Test;
public class SymmLQTest {
public void saundersTest(final int n, final boolean goodb,
final boolean precon, final double shift,
final double pertbn) {
final RealLinearOperator a = new RealLinearOperator() {
@Override
public RealVector operate(final RealVector x) {
if (x.getDimension() != n) {
throw new DimensionMismatchException(x.getDimension(), n);
}
final double[] y = new double[n];
for (int i = 0; i < n; i++) {
y[i] = (i + 1) * 1.1 / n * x.getEntry(i);
}
return new ArrayRealVector(y, false);
}
@Override
public int getRowDimension() {
return n;
}
@Override
public int getColumnDimension() {
return n;
}
};
final double shiftm = shift;
final double pertm = FastMath.abs(pertbn);
final RealLinearOperator minv;
if (precon) {
minv = new RealLinearOperator() {
@Override
public int getRowDimension() {
return n;
}
@Override
public int getColumnDimension() {
return n;
}
@Override
public RealVector operate(final RealVector x) {
if (x.getDimension() != n) {
throw new DimensionMismatchException(x.getDimension(),
n);
}
final double[] y = new double[n];
for (int i = 0; i < n; i++) {
double d = (i + 1) * 1.1 / n;
d = FastMath.abs(d - shiftm);
if (i % 10 == 0) {
d += pertm;
}
y[i] = x.getEntry(i) / d;
}
return new ArrayRealVector(y, false);
}
};
} else {
minv = null;
}
final RealVector xtrue = new ArrayRealVector(n);
for (int i = 0; i < n; i++) {
xtrue.setEntry(i, n - i);
}
final RealVector b = a.operate(xtrue);
b.combineToSelf(1.0, -shift, xtrue);
final SymmLQ solver = new SymmLQ(2 * n, 1E-12, true);
final RealVector x = solver.solve(a, minv, b, goodb, shift);
final RealVector y = a.operate(x);
final RealVector r1 = new ArrayRealVector(n);
for (int i = 0; i < n; i++) {
final double bi = b.getEntry(i);
final double yi = y.getEntry(i);
final double xi = x.getEntry(i);
r1.setEntry(i, bi - yi + shift * xi);
}
final double enorm = x.subtract(xtrue).getNorm() / xtrue.getNorm();
final double etol = 1E-5;
Assert.assertTrue("enorm=" + enorm + ", " +
solver.getIterationManager().getIterations(), enorm <= etol);
}
@Test
public void testSolveSaunders1() {
saundersTest(1, false, false, 0., 0.);
}
@Test
public void testSolveSaunders2() {
saundersTest(2, false, false, 0., 0.);
}
@Test
public void testSolveSaunders3() {
saundersTest(1, false, true, 0., 0.);
}
@Test
public void testSolveSaunders4() {
saundersTest(2, false, true, 0., 0.);
}
@Test
public void testSolveSaunders5() {
saundersTest(5, false, true, 0., 0.);
}
@Test
public void testSolveSaunders6() {
saundersTest(5, false, true, 0.25, 0.);
}
@Test
public void testSolveSaunders7() {
saundersTest(50, false, false, 0., 0.);
}
@Test
public void testSolveSaunders8() {
saundersTest(50, false, false, 0.25, 0.);
}
@Test
public void testSolveSaunders9() {
saundersTest(50, false, true, 0., 0.10);
}
@Test
public void testSolveSaunders10() {
saundersTest(50, false, true, 0.25, 0.10);
}
@Test
public void testSolveSaunders11() {
saundersTest(1, true, false, 0., 0.);
}
@Test
public void testSolveSaunders12() {
saundersTest(2, true, false, 0., 0.);
}
@Test
public void testSolveSaunders13() {
saundersTest(1, true, true, 0., 0.);
}
@Test
public void testSolveSaunders14() {
saundersTest(2, true, true, 0., 0.);
}
@Test
public void testSolveSaunders15() {
saundersTest(5, true, true, 0., 0.);
}
@Test
public void testSolveSaunders16() {
saundersTest(5, true, true, 0.25, 0.);
}
@Test
public void testSolveSaunders17() {
saundersTest(50, true, false, 0., 0.);
}
@Test
public void testSolveSaunders18() {
saundersTest(50, true, false, 0.25, 0.);
}
@Test
public void testSolveSaunders19() {
saundersTest(50, true, true, 0., 0.10);
}
@Test
public void testSolveSaunders20() {
saundersTest(50, true, true, 0.25, 0.10);
}
@Test(expected = NonSquareOperatorException.class)
public void testNonSquareOperator() {
final Array2DRowRealMatrix a = new Array2DRowRealMatrix(2, 3);
final IterativeLinearSolver solver;
solver = new SymmLQ(10, 0., false);
final ArrayRealVector b = new ArrayRealVector(a.getRowDimension());
final ArrayRealVector x = new ArrayRealVector(a.getColumnDimension());
solver.solve(a, b, x);
}
@Test(expected = DimensionMismatchException.class)
public void testDimensionMismatchRightHandSide() {
final Array2DRowRealMatrix a = new Array2DRowRealMatrix(3, 3);
final IterativeLinearSolver solver;
solver = new SymmLQ(10, 0., false);
final ArrayRealVector b = new ArrayRealVector(2);
solver.solve(a, b);
}
@Test(expected = DimensionMismatchException.class)
public void testDimensionMismatchSolution() {
final Array2DRowRealMatrix a = new Array2DRowRealMatrix(3, 3);
final IterativeLinearSolver solver;
solver = new SymmLQ(10, 0., false);
final ArrayRealVector b = new ArrayRealVector(3);
final ArrayRealVector x = new ArrayRealVector(2);
solver.solve(a, b, x);
}
@Test
public void testUnpreconditionedSolution() {
final int n = 5;
final int maxIterations = 100;
final RealLinearOperator a = new HilbertMatrix(n);
final InverseHilbertMatrix ainv = new InverseHilbertMatrix(n);
final IterativeLinearSolver solver;
solver = new SymmLQ(maxIterations, 1E-10, true);
final RealVector b = new ArrayRealVector(n);
for (int j = 0; j < n; j++) {
b.set(0.);
b.setEntry(j, 1.);
final RealVector x = solver.solve(a, b);
for (int i = 0; i < n; i++) {
final double actual = x.getEntry(i);
final double expected = ainv.getEntry(i, j);
final double delta = 1E-6 * FastMath.abs(expected);
final String msg = String.format("entry[%d][%d]", i, j);
Assert.assertEquals(msg, expected, actual, delta);
}
}
}
@Test
public void testUnpreconditionedInPlaceSolutionWithInitialGuess() {
final int n = 5;
final int maxIterations = 100;
final RealLinearOperator a = new HilbertMatrix(n);
final InverseHilbertMatrix ainv = new InverseHilbertMatrix(n);
final IterativeLinearSolver solver;
solver = new SymmLQ(maxIterations, 1E-10, true);
final RealVector b = new ArrayRealVector(n);
for (int j = 0; j < n; j++) {
b.set(0.);
b.setEntry(j, 1.);
final RealVector x0 = new ArrayRealVector(n);
x0.set(1.);
final RealVector x = solver.solveInPlace(a, b, x0);
Assert.assertSame("x should be a reference to x0", x0, x);
for (int i = 0; i < n; i++) {
final double actual = x.getEntry(i);
final double expected = ainv.getEntry(i, j);
final double delta = 1E-6 * FastMath.abs(expected);
final String msg = String.format("entry[%d][%d)", i, j);
Assert.assertEquals(msg, expected, actual, delta);
}
}
}
@Test
public void testUnpreconditionedSolutionWithInitialGuess() {
final int n = 5;
final int maxIterations = 100;
final RealLinearOperator a = new HilbertMatrix(n);
final InverseHilbertMatrix ainv = new InverseHilbertMatrix(n);
final IterativeLinearSolver solver;
solver = new SymmLQ(maxIterations, 1E-10, true);
final RealVector b = new ArrayRealVector(n);
for (int j = 0; j < n; j++) {
b.set(0.);
b.setEntry(j, 1.);
final RealVector x0 = new ArrayRealVector(n);
x0.set(1.);
final RealVector x = solver.solve(a, b, x0);
Assert.assertNotSame("x should not be a reference to x0", x0, x);
for (int i = 0; i < n; i++) {
final double actual = x.getEntry(i);
final double expected = ainv.getEntry(i, j);
final double delta = 1E-6 * FastMath.abs(expected);
final String msg = String.format("entry[%d][%d]", i, j);
Assert.assertEquals(msg, expected, actual, delta);
Assert.assertEquals(msg, x0.getEntry(i), 1., Math.ulp(1.));
}
}
}
@Test(expected = NonSquareOperatorException.class)
public void testNonSquarePreconditioner() {
final Array2DRowRealMatrix a = new Array2DRowRealMatrix(2, 2);
final RealLinearOperator m = new RealLinearOperator() {
@Override
public RealVector operate(final RealVector x) {
throw new UnsupportedOperationException();
}
@Override
public int getRowDimension() {
return 2;
}
@Override
public int getColumnDimension() {
return 3;
}
};
final PreconditionedIterativeLinearSolver solver;
solver = new SymmLQ(10, 0., false);
final ArrayRealVector b = new ArrayRealVector(a.getRowDimension());
solver.solve(a, m, b);
}
@Test(expected = DimensionMismatchException.class)
public void testMismatchedOperatorDimensions() {
final Array2DRowRealMatrix a = new Array2DRowRealMatrix(2, 2);
final RealLinearOperator m = new RealLinearOperator() {
@Override
public RealVector operate(final RealVector x) {
throw new UnsupportedOperationException();
}
@Override
public int getRowDimension() {
return 3;
}
@Override
public int getColumnDimension() {
return 3;
}
};
final PreconditionedIterativeLinearSolver solver;
solver = new SymmLQ(10, 0d, false);
final ArrayRealVector b = new ArrayRealVector(a.getRowDimension());
solver.solve(a, m, b);
}
@Test(expected = NonPositiveDefiniteOperatorException.class)
public void testNonPositiveDefinitePreconditioner() {
final Array2DRowRealMatrix a = new Array2DRowRealMatrix(2, 2);
a.setEntry(0, 0, 1d);
a.setEntry(0, 1, 2d);
a.setEntry(1, 0, 3d);
a.setEntry(1, 1, 4d);
final RealLinearOperator m = new RealLinearOperator() {
@Override
public RealVector operate(final RealVector x) {
final ArrayRealVector y = new ArrayRealVector(2);
y.setEntry(0, -x.getEntry(0));
y.setEntry(1, -x.getEntry(1));
return y;
}
@Override
public int getRowDimension() {
return 2;
}
@Override
public int getColumnDimension() {
return 2;
}
};
final PreconditionedIterativeLinearSolver solver;
solver = new SymmLQ(10, 0d, true);
final ArrayRealVector b = new ArrayRealVector(2);
b.setEntry(0, -1d);
b.setEntry(1, -1d);
solver.solve(a, m, b);
}
@Test
public void testPreconditionedSolution() {
final int n = 8;
final int maxIterations = 100;
final RealLinearOperator a = new HilbertMatrix(n);
final InverseHilbertMatrix ainv = new InverseHilbertMatrix(n);
final RealLinearOperator m = JacobiPreconditioner.create(a);
final PreconditionedIterativeLinearSolver solver;
solver = new SymmLQ(maxIterations, 1E-15, true);
final RealVector b = new ArrayRealVector(n);
for (int j = 0; j < n; j++) {
b.set(0.);
b.setEntry(j, 1.);
final RealVector x = solver.solve(a, m, b);
for (int i = 0; i < n; i++) {
final double actual = x.getEntry(i);
final double expected = ainv.getEntry(i, j);
final double delta = 1E-6 * FastMath.abs(expected);
final String msg = String.format("coefficient (%d, %d)", i, j);
Assert.assertEquals(msg, expected, actual, delta);
}
}
}
@Test
public void testPreconditionedSolution2() {
final int n = 100;
final int maxIterations = 100000;
final Array2DRowRealMatrix a = new Array2DRowRealMatrix(n, n);
double daux = 1.;
for (int i = 0; i < n; i++) {
a.setEntry(i, i, daux);
daux *= 1.2;
for (int j = i + 1; j < n; j++) {
if (i == j) {
} else {
final double value = 1.0;
a.setEntry(i, j, value);
a.setEntry(j, i, value);
}
}
}
final RealLinearOperator m = JacobiPreconditioner.create(a);
final PreconditionedIterativeLinearSolver prec;
final IterativeLinearSolver unprec;
prec = new SymmLQ(maxIterations, 1E-15, true);
unprec = new SymmLQ(maxIterations, 1E-15, true);
final RealVector b = new ArrayRealVector(n);
final String pattern = "preconditioned SymmLQ (%d iterations) should"
+ " have been faster than unpreconditioned (%d iterations)";
String msg;
for (int j = 0; j < 1; j++) {
b.set(0.);
b.setEntry(j, 1.);
final RealVector px = prec.solve(a, m, b);
final RealVector x = unprec.solve(a, b);
final int np = prec.getIterationManager().getIterations();
final int nup = unprec.getIterationManager().getIterations();
msg = String.format(pattern, np, nup);
for (int i = 0; i < n; i++) {
msg = String.format("row %d, column %d", i, j);
final double expected = x.getEntry(i);
final double actual = px.getEntry(i);
final double delta = 5E-5 * FastMath.abs(expected);
Assert.assertEquals(msg, expected, actual, delta);
}
}
}
@Test
public void testEventManagement() {
final int n = 5;
final int maxIterations = 100;
final RealLinearOperator a = new HilbertMatrix(n);
final IterativeLinearSolver solver;
/*
* count[0] = number of calls to initializationPerformed
* count[1] = number of calls to iterationStarted
* count[2] = number of calls to iterationPerformed
* count[3] = number of calls to terminationPerformed
*/
final int[] count = new int[] {0, 0, 0, 0};
final RealVector xFromListener = new ArrayRealVector(n);
final IterationListener listener = new IterationListener() {
@Override
public void initializationPerformed(final IterationEvent e) {
++count[0];
}
@Override
public void iterationPerformed(final IterationEvent e) {
++count[2];
Assert.assertEquals("iteration performed",
count[2],
e.getIterations() - 1);
}
@Override
public void iterationStarted(final IterationEvent e) {
++count[1];
Assert.assertEquals("iteration started",
count[1],
e.getIterations() - 1);
}
@Override
public void terminationPerformed(final IterationEvent e) {
++count[3];
final IterativeLinearSolverEvent ilse;
ilse = (IterativeLinearSolverEvent) e;
xFromListener.setSubVector(0, ilse.getSolution());
}
};
solver = new SymmLQ(maxIterations, 1E-10, true);
solver.getIterationManager().addIterationListener(listener);
final RealVector b = new ArrayRealVector(n);
for (int j = 0; j < n; j++) {
Arrays.fill(count, 0);
b.set(0.);
b.setEntry(j, 1.);
final RealVector xFromSolver = solver.solve(a, b);
String msg = String.format("column %d (initialization)", j);
Assert.assertEquals(msg, 1, count[0]);
msg = String.format("column %d (finalization)", j);
Assert.assertEquals(msg, 1, count[3]);
/*
* Check that solution is not "over-refined". When the last
* iteration has occurred, no further refinement should be
* performed.
*/
for (int i = 0; i < n; i++){
msg = String.format("row %d, column %d", i, j);
final double expected = xFromSolver.getEntry(i);
final double actual = xFromListener.getEntry(i);
Assert.assertEquals(msg, expected, actual, 0.0);
}
}
}
@Test(expected = NonSelfAdjointOperatorException.class)
public void testNonSelfAdjointOperator() {
final RealLinearOperator a;
a = new Array2DRowRealMatrix(new double[][] {
{1., 2., 3.},
{2., 4., 5.},
{2.999, 5., 6.}
});
final RealVector b;
b = new ArrayRealVector(new double[] {1., 1., 1.});
new SymmLQ(100, 1., true).solve(a, b);
}
@Test(expected = NonSelfAdjointOperatorException.class)
public void testNonSelfAdjointPreconditioner() {
final RealLinearOperator a = new Array2DRowRealMatrix(new double[][] {
{1., 2., 3.},
{2., 4., 5.},
{3., 5., 6.}
});
final Array2DRowRealMatrix mMat;
mMat = new Array2DRowRealMatrix(new double[][] {
{1., 0., 1.},
{0., 1., 0.},
{0., 0., 1.}
});
final DecompositionSolver mSolver;
mSolver = new LUDecomposition(mMat).getSolver();
final RealLinearOperator minv = new RealLinearOperator() {
@Override
public RealVector operate(final RealVector x) {
return mSolver.solve(x);
}
@Override
public int getRowDimension() {
return mMat.getRowDimension();
}
@Override
public int getColumnDimension() {
return mMat.getColumnDimension();
}
};
final RealVector b = new ArrayRealVector(new double[] {
1., 1., 1.
});
new SymmLQ(100, 1., true).solve(a, minv, b);
}
@Test
public void testUnpreconditionedNormOfResidual() {
final int n = 5;
final int maxIterations = 100;
final RealLinearOperator a = new HilbertMatrix(n);
final IterativeLinearSolver solver;
final IterationListener listener = new IterationListener() {
private void doTestNormOfResidual(final IterationEvent e) {
final IterativeLinearSolverEvent evt;
evt = (IterativeLinearSolverEvent) e;
final RealVector x = evt.getSolution();
final RealVector b = evt.getRightHandSideVector();
final RealVector r = b.subtract(a.operate(x));
final double rnorm = r.getNorm();
Assert.assertEquals("iteration performed (residual)",
rnorm, evt.getNormOfResidual(),
FastMath.max(1E-5 * rnorm, 1E-10));
}
@Override
public void initializationPerformed(final IterationEvent e) {
doTestNormOfResidual(e);
}
@Override
public void iterationPerformed(final IterationEvent e) {
doTestNormOfResidual(e);
}
@Override
public void iterationStarted(final IterationEvent e) {
doTestNormOfResidual(e);
}
@Override
public void terminationPerformed(final IterationEvent e) {
doTestNormOfResidual(e);
}
};
solver = new SymmLQ(maxIterations, 1E-10, true);
solver.getIterationManager().addIterationListener(listener);
final RealVector b = new ArrayRealVector(n);
for (int j = 0; j < n; j++) {
b.set(0.);
b.setEntry(j, 1.);
solver.solve(a, b);
}
}
@Test
public void testPreconditionedNormOfResidual() {
final int n = 5;
final int maxIterations = 100;
final RealLinearOperator a = new HilbertMatrix(n);
final JacobiPreconditioner m = JacobiPreconditioner.create(a);
final RealLinearOperator p = m.sqrt();
final PreconditionedIterativeLinearSolver solver;
final IterationListener listener = new IterationListener() {
private void doTestNormOfResidual(final IterationEvent e) {
final IterativeLinearSolverEvent evt;
evt = (IterativeLinearSolverEvent) e;
final RealVector x = evt.getSolution();
final RealVector b = evt.getRightHandSideVector();
final RealVector r = b.subtract(a.operate(x));
final double rnorm = p.operate(r).getNorm();
Assert.assertEquals("iteration performed (residual)",
rnorm, evt.getNormOfResidual(),
FastMath.max(1E-5 * rnorm, 1E-10));
}
@Override
public void initializationPerformed(final IterationEvent e) {
doTestNormOfResidual(e);
}
@Override
public void iterationPerformed(final IterationEvent e) {
doTestNormOfResidual(e);
}
@Override
public void iterationStarted(final IterationEvent e) {
doTestNormOfResidual(e);
}
@Override
public void terminationPerformed(final IterationEvent e) {
doTestNormOfResidual(e);
}
};
solver = new SymmLQ(maxIterations, 1E-10, true);
solver.getIterationManager().addIterationListener(listener);
final RealVector b = new ArrayRealVector(n);
for (int j = 0; j < n; j++) {
b.set(0.);
b.setEntry(j, 1.);
solver.solve(a, m, b);
}
}
}