/*
* 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.mahout.math;
import java.util.List;
import com.google.common.collect.Lists;
import org.apache.mahout.math.function.DoubleDoubleFunction;
import org.apache.mahout.math.function.Functions;
import org.apache.mahout.math.stats.OnlineSummarizer;
import org.junit.Ignore;
import org.junit.Test;
public final class QRDecompositionTest extends MahoutTestCase {
@Test
public void randomMatrix() {
Matrix a = new DenseMatrix(60, 60).assign(Functions.random());
QRDecomposition qr = new QRDecomposition(a);
// how close is Q to actually being orthornormal?
double maxIdent = qr.getQ().transpose().times(qr.getQ()).viewDiagonal().assign(Functions.plus(-1)).norm(1);
assertEquals(0, maxIdent, 1.0e-13);
// how close is Q R to the original value of A?
Matrix z = qr.getQ().times(qr.getR()).minus(a);
double maxError = z.aggregate(Functions.MIN, Functions.ABS);
assertEquals(0, maxError, 1.0e-13);
}
@Test
public void rank1() {
Matrix x = new DenseMatrix(3, 3);
x.viewRow(0).assign(new double[]{1, 2, 3});
x.viewRow(1).assign(new double[]{2, 4, 6});
x.viewRow(2).assign(new double[]{3, 6, 9});
QRDecomposition qr = new QRDecomposition(x);
assertFalse(qr.hasFullRank());
assertEquals(0, new DenseVector(new double[]{3.741657, 7.483315, 11.22497}).aggregate(qr.getR().viewRow(0), Functions.PLUS, new DoubleDoubleFunction() {
@Override
public double apply(double arg1, double arg2) {
return Math.abs(arg1) - Math.abs(arg2);
}
}), 1.0e-5);
}
@Test
public void fullRankTall() {
Matrix x = matrix();
QRDecomposition qr = new QRDecomposition(x);
assertTrue(qr.hasFullRank());
Matrix rRef = reshape(new double[]{
-2.99129686445138, 0, 0, 0, 0,
-0.0282260628674372, -2.38850244769059, 0, 0, 0,
0.733739310355871, 1.48042000631646, 2.29051263117895, 0, 0,
-0.0394082168269326, 0.282829484207801, -0.00438521041803086, -2.90823198084203, 0,
0.923669647838536, 1.76679276072492, 0.637690104222683, -0.225890909498753, -1.35732293800944},
5, 5);
Matrix r = qr.getR();
// check identity down to sign
assertEquals(0, r.clone().assign(Functions.ABS).minus(rRef.clone().assign(Functions.ABS)).aggregate(Functions.PLUS, Functions.IDENTITY), 1.0e-12);
Matrix qRef = reshape(new double[]{
-0.165178287646573, 0.0510035857637869, 0.13985915987379, -0.120173729496501,
-0.453198314345324, 0.644400679630493, -0.503117990820608, 0.24968739845381,
0.323968339146224, -0.465266080134262, 0.276508948773268, -0.687909700644343,
0.0544048888907195, -0.0166677718378263, 0.171309755790717, 0.310339001630029,
0.674790532821663, 0.0058166082200493, -0.381707516461884, 0.300504956413142,
-0.105751091334003, 0.410450870871096, 0.31113446615821, 0.179338172684956,
0.361951807617901, 0.763921725548796, 0.380327892605634, -0.287274944594054,
0.0311604042556675, 0.0386096858143961, 0.0387156960650472, -0.232975755728917,
0.0358178276684149, 0.173105775703199, 0.327321867815603, 0.328671945345279,
-0.36015879836344, -0.444261660176044, 0.09438499563253, 0.646216148583769
}, 8, 5);
printMatrix("qRef", qRef);
Matrix q = qr.getQ();
printMatrix("q", q);
assertEquals(0, q.clone().assign(Functions.ABS).minus(qRef.clone().assign(Functions.ABS)).aggregate(Functions.PLUS, Functions.IDENTITY), 1.0e-12);
Matrix x1 = qr.solve(reshape(new double[]{
-0.0178247686747641, 0.68631714634098, -0.335464858468858, 1.50249941751569,
-0.669901640772149, -0.977025038942455, -1.18857546169856, -1.24792900492054
}, 8, 1));
Matrix xref = reshape(new double[]{
-0.0127440093664874, 0.655825940180799, -0.100755415991702, -0.0349559562697406,
-0.190744297762028
}, 5, 1);
printMatrix("x1", x1);
printMatrix("xref", xref);
assertEquals(xref, x1, 1.0e-8);
}
@Test
public void fullRankWide() {
Matrix x = matrix().transpose();
QRDecomposition qr = new QRDecomposition(x);
assertTrue(qr.hasFullRank());
Matrix rActual = qr.getR();
Matrix rRef = reshape(new double[]{
-2.42812464965842, 0, 0, 0, 0,
0.303587286111356, -2.91663643494775, 0, 0, 0,
-0.201812474153156, -0.765485720168378, 1.09989373598954, 0, 0,
1.47980701097885, -0.637545820524326, -1.55519859337935, 0.844655127991726, 0,
0.0248883129453161, 0.00115010570270549, -0.236340588891252, -0.092924118200147, 1.42910099545547,
-1.1678472412429, 0.531245845248056, 0.351978196071514, -1.03241474816555, -2.20223861735426,
-0.887809959067632, 0.189731251982918, -0.504321849233586, 0.490484123999836, 1.21266692336743,
-0.633888169775463, 1.04738559065986, 0.284041239547031, 0.578183510077156, -0.942314870832456
}, 5, 8);
printMatrix("rRef", rRef);
printMatrix("rActual", rActual);
assertEquals(0, rActual.clone().assign(Functions.ABS).minus(rRef.clone().assign(Functions.ABS)).aggregate(Functions.PLUS, Functions.IDENTITY), 1.0e-12);
// assertEquals(rRef, rActual, 1.0e-8);
Matrix qRef = reshape(new double[]{
-0.203489262374627, 0.316761677948356, -0.784155643293468, 0.394321494579, -0.29641971170211,
0.0311283614803723, -0.34755265020736, 0.137138511478328, 0.848579887681972, 0.373287266507375,
-0.39603700561249, -0.787812566647329, -0.377864833067864, -0.275080943427399, 0.0636764674878229,
0.0763976893309043, -0.318551137554327, 0.286407036668598, 0.206004127289883, -0.876482672226889,
0.89159476695423, -0.238213616975551, -0.376141107880836, -0.0794701657055114, 0.0227025098210165
}, 5, 5);
Matrix q = qr.getQ();
printMatrix("qRef", qRef);
printMatrix("q", q);
assertEquals(0, q.clone().assign(Functions.ABS).minus(qRef.clone().assign(Functions.ABS)).aggregate(Functions.PLUS, Functions.IDENTITY), 1.0e-12);
// assertEquals(qRef, q, 1.0e-8);
Matrix x1 = qr.solve(b());
Matrix xRef = reshape(new double[]{
-0.182580239668147, -0.437233627652114, 0.138787653097464, 0.672934739896228, -0.131420217069083, 0, 0, 0
}, 8, 1);
printMatrix("xRef", xRef);
printMatrix("x", x1);
assertEquals(xRef, x1, 1.0e-8);
assertEquals(x, qr.getQ().times(qr.getR()), 1.0e-15);
}
// TODO: the speedup constant should be checked and oddly, the times don't increase as the counts increase
@Ignore
public void fasterThanBefore() {
OnlineSummarizer s1 = new OnlineSummarizer();
OnlineSummarizer s2 = new OnlineSummarizer();
Matrix a = new DenseMatrix(60, 60).assign(Functions.random());
decompositionSpeedCheck(new Decomposer() {
@Override
public QR decompose(Matrix a) {
return new QRDecomposition(a);
}
}, s1, a, "new");
decompositionSpeedCheck(new Decomposer() {
@Override
public QR decompose(Matrix a) {
return new OldQRDecomposition(a);
}
}, s2, a, "old");
// should be much more than twice as fast.
System.out.printf("Speedup is about %.1f times\n", s2.getMedian() / s1.getMedian());
assertTrue(s1.getMedian() < 0.5 * s2.getMedian());
}
private interface Decomposer {
QR decompose(Matrix a);
}
private static void decompositionSpeedCheck(Decomposer qrf, OnlineSummarizer s1, Matrix a, String label) {
int n = 0;
List<Integer> counts = Lists.newArrayList(10, 20, 50, 100, 200, 500);
for (int k : counts) {
double warmup = 0;
double other = 0;
n += k;
for (int i = 0; i < k; i++) {
QR qr = qrf.decompose(a);
warmup = Math.max(warmup, qr.getQ().transpose().times(qr.getQ()).viewDiagonal().assign(Functions.plus(-1)).norm(1));
Matrix z = qr.getQ().times(qr.getR()).minus(a);
other = Math.max(other, z.aggregate(Functions.MIN, Functions.ABS));
}
double maxIdent = 0;
double maxError = 0;
long t0 = System.nanoTime();
for (int i = 0; i < n; i++) {
QR qr = qrf.decompose(a);
maxIdent = Math.max(maxIdent, qr.getQ().transpose().times(qr.getQ()).viewDiagonal().assign(Functions.plus(-1)).norm(1));
Matrix z = qr.getQ().times(qr.getR()).minus(a);
maxError = Math.max(maxError, z.aggregate(Functions.MIN, Functions.ABS));
}
long t1 = System.nanoTime();
if (k > 100) {
s1.add(t1 - t0);
}
System.out.printf("%s %d\t%.1f\t%g\t%g\t%g\n", label, n, (t1 - t0) / 1.0e3 / n, maxIdent, maxError, warmup);
}
}
private static void assertEquals(Matrix ref, Matrix actual, double epsilon) {
assertEquals(0, ref.minus(actual).aggregate(Functions.MAX, Functions.ABS), epsilon);
}
private static void printMatrix(String name, Matrix m) {
int rows = m.numRows();
int columns = m.numCols();
System.out.printf("%s - %d x %d\n", name, rows, columns);
for (int i = 0; i < rows; i++) {
for (int j = 0; j < columns; j++) {
System.out.printf("%10.5f", m.get(i, j));
}
System.out.printf("\n");
}
System.out.printf("\n");
System.out.printf("\n");
}
private static Matrix matrix() {
double[] values = {
0.494097293912641, -0.152566866170993, -0.418360266395271, 0.359475300232312,
1.35565069667582, -1.92759373242903, 1.50497526839076, -0.746889132087904,
-0.769136838293565, 1.10984954080986, -0.664389974392489, 1.6464660350229,
-0.11715420616969, 0.0216221197371269, -0.394972730980765, -0.748293157213142,
1.90402764664962, -0.638042862848559, -0.362336344669668, -0.418261074380526,
-0.494211543128429, 1.38828971158414, 0.597110366867923, 1.05341387608687,
-0.957461740877418, -2.35528802598249, -1.03171458944128, 0.644319090271635,
-0.0569108993041965, -0.14419465550881, -0.0456801828174936,
0.754694392571835, 0.719744008628535, -1.17873249802301, -0.155887528905918,
-1.5159868405466, 0.0918931582603128, 1.42179027361583, -0.100495054250176,
0.0687986548485584
};
return reshape(values, 8, 5);
}
private static Matrix reshape(double[] values, int rows, int columns) {
Matrix m = new DenseMatrix(rows, columns);
int i = 0;
for (double v : values) {
m.set(i % rows, i / rows, v);
i++;
}
return m;
}
private static Matrix b() {
return reshape(new double[]
{-0.0178247686747641, 0.68631714634098, -0.335464858468858, 1.50249941751569, -0.669901640772149}, 5, 1);
}
}