summaryrefslogtreecommitdiff
path: root/src/de/lmu/ifi/dbs/elki/math/linearalgebra/Matrix.java
diff options
context:
space:
mode:
Diffstat (limited to 'src/de/lmu/ifi/dbs/elki/math/linearalgebra/Matrix.java')
-rw-r--r--src/de/lmu/ifi/dbs/elki/math/linearalgebra/Matrix.java526
1 files changed, 239 insertions, 287 deletions
diff --git a/src/de/lmu/ifi/dbs/elki/math/linearalgebra/Matrix.java b/src/de/lmu/ifi/dbs/elki/math/linearalgebra/Matrix.java
index 3f513720..6826f67d 100644
--- a/src/de/lmu/ifi/dbs/elki/math/linearalgebra/Matrix.java
+++ b/src/de/lmu/ifi/dbs/elki/math/linearalgebra/Matrix.java
@@ -4,7 +4,7 @@ package de.lmu.ifi.dbs.elki.math.linearalgebra;
This file is part of ELKI:
Environment for Developing KDD-Applications Supported by Index-Structures
- Copyright (C) 2013
+ Copyright (C) 2014
Ludwig-Maximilians-Universität München
Lehr- und Forschungseinheit für Datenbanksysteme
ELKI Development Team
@@ -31,7 +31,6 @@ import java.util.ArrayList;
import java.util.Arrays;
import java.util.logging.Logger;
-import de.lmu.ifi.dbs.elki.data.RationalNumber;
import de.lmu.ifi.dbs.elki.logging.LoggingConfiguration;
import de.lmu.ifi.dbs.elki.math.MathUtil;
import de.lmu.ifi.dbs.elki.utilities.FormatUtil;
@@ -58,9 +57,15 @@ public class Matrix {
public static final double DELTA = 1E-3;
/**
- * Error: matrix not square.
+ * Small value to increment diagonally of a matrix in order to avoid
+ * singularity before building the inverse.
*/
- public static final String ERR_NOTSQUARE = "All rows must have the same length.";
+ public static final double SINGULARITY_CHEAT = 1E-9;
+
+ /**
+ * Error: matrix not rectangular.
+ */
+ public static final String ERR_NOTRECTANGULAR = "All rows must have the same length.";
/**
* Error: matrix indexes incorrect
@@ -112,8 +117,8 @@ public class Matrix {
public Matrix(final int m, final int n, final double s) {
this.columndimension = n;
elements = new double[m][n];
- for (int i = 0; i < m; i++) {
- for (int j = 0; j < n; j++) {
+ for(int i = 0; i < m; i++) {
+ for(int j = 0; j < n; j++) {
elements[i][j] = s;
}
}
@@ -128,31 +133,15 @@ public class Matrix {
*/
public Matrix(final double[][] elements) {
columndimension = elements[0].length;
- for (int i = 0; i < elements.length; i++) {
- if (elements[i].length != columndimension) {
- throw new IllegalArgumentException(ERR_NOTSQUARE);
+ for(int i = 0; i < elements.length; i++) {
+ if(elements[i].length != columndimension) {
+ throw new IllegalArgumentException(ERR_NOTRECTANGULAR);
}
}
this.elements = elements;
}
/**
- * Constructs a Matrix for a given array of arrays of {@link RationalNumber}s.
- *
- * @param q an array of arrays of RationalNumbers. q is not checked for
- * consistency (i.e. whether all rows are of equal length)
- */
- public Matrix(final RationalNumber[][] q) {
- columndimension = q[0].length;
- elements = new double[q.length][columndimension];
- for (int row = 0; row < q.length; row++) {
- for (int col = 0; col < q[row].length; col++) {
- elements[row][col] = q[row][col].doubleValue();
- }
- }
- }
-
- /**
* Construct a matrix from a one-dimensional packed array
*
* @param values One-dimensional array of doubles, packed by columns (ala
@@ -162,12 +151,12 @@ public class Matrix {
*/
public Matrix(final double values[], final int m) {
columndimension = (m != 0 ? values.length / m : 0);
- if (m * columndimension != values.length) {
+ if(m * columndimension != values.length) {
throw new IllegalArgumentException("Array length must be a multiple of m.");
}
elements = new double[m][columndimension];
- for (int i = 0; i < m; i++) {
- for (int j = 0; j < columndimension; j++) {
+ for(int i = 0; i < m; i++) {
+ for(int j = 0; j < columndimension; j++) {
elements[i][j] = values[i + j * m];
}
}
@@ -193,9 +182,9 @@ public class Matrix {
final int m = A.length;
final int n = A[0].length;
final Matrix X = new Matrix(m, n);
- for (int i = 0; i < m; i++) {
- if (A[i].length != n) {
- throw new IllegalArgumentException(ERR_NOTSQUARE);
+ for(int i = 0; i < m; i++) {
+ if(A[i].length != n) {
+ throw new IllegalArgumentException(ERR_NOTRECTANGULAR);
}
System.arraycopy(A[i], 0, X.elements[i], 0, n);
}
@@ -210,7 +199,7 @@ public class Matrix {
*/
public static final Matrix unitMatrix(final int dim) {
final double[][] e = new double[dim][dim];
- for (int i = 0; i < dim; i++) {
+ for(int i = 0; i < dim; i++) {
e[i][i] = 1;
}
return new Matrix(e);
@@ -236,8 +225,8 @@ public class Matrix {
*/
public static final Matrix random(final int m, final int n) {
final Matrix A = new Matrix(m, n);
- for (int i = 0; i < m; i++) {
- for (int j = 0; j < n; j++) {
+ for(int i = 0; i < m; i++) {
+ for(int j = 0; j < n; j++) {
A.elements[i][j] = Math.random();
}
}
@@ -253,7 +242,7 @@ public class Matrix {
*/
public static final Matrix identity(final int m, final int n) {
final Matrix A = new Matrix(m, n);
- for (int i = 0; i < Math.min(m, n); i++) {
+ for(int i = 0; i < Math.min(m, n); i++) {
A.elements[i][i] = 1.0;
}
return A;
@@ -268,7 +257,7 @@ public class Matrix {
*/
public static final Matrix diagonal(final double[] diagonal) {
final Matrix result = new Matrix(diagonal.length, diagonal.length);
- for (int i = 0; i < diagonal.length; i++) {
+ for(int i = 0; i < diagonal.length; i++) {
result.elements[i][i] = diagonal[i];
}
return result;
@@ -283,7 +272,7 @@ public class Matrix {
*/
public static final Matrix diagonal(final Vector diagonal) {
final Matrix result = new Matrix(diagonal.elements.length, diagonal.elements.length);
- for (int i = 0; i < diagonal.elements.length; i++) {
+ for(int i = 0; i < diagonal.elements.length; i++) {
result.elements[i][i] = diagonal.elements[i];
}
return result;
@@ -296,7 +285,7 @@ public class Matrix {
*/
public final Matrix copy() {
final Matrix X = new Matrix(elements.length, columndimension);
- for (int i = 0; i < elements.length; i++) {
+ for(int i = 0; i < elements.length; i++) {
System.arraycopy(elements[i], 0, X.elements[i], 0, columndimension);
}
return X;
@@ -326,7 +315,7 @@ public class Matrix {
*/
public final double[][] getArrayCopy() {
final double[][] C = new double[elements.length][];
- for (int i = 0; i < elements.length; i++) {
+ for(int i = 0; i < elements.length; i++) {
C[i] = elements[i].clone();
}
return C;
@@ -397,7 +386,7 @@ public class Matrix {
*/
public final double[] getRowPackedCopy() {
double[] vals = new double[elements.length * columndimension];
- for (int i = 0; i < elements.length; i++) {
+ for(int i = 0; i < elements.length; i++) {
System.arraycopy(elements[i], 0, vals, i * columndimension, columndimension);
}
return vals;
@@ -410,8 +399,8 @@ public class Matrix {
*/
public final double[] getColumnPackedCopy() {
final double[] vals = new double[elements.length * columndimension];
- for (int i = 0; i < elements.length; i++) {
- for (int j = 0; j < columndimension; j++) {
+ for(int i = 0; i < elements.length; i++) {
+ for(int j = 0; j < columndimension; j++) {
vals[i + j * elements.length] = elements[i][j];
}
}
@@ -431,10 +420,11 @@ public class Matrix {
public final Matrix getMatrix(final int i0, final int i1, final int j0, final int j1) {
final Matrix X = new Matrix(i1 - i0 + 1, j1 - j0 + 1);
try {
- for (int i = i0; i <= i1; i++) {
+ for(int i = i0; i <= i1; i++) {
System.arraycopy(elements[i], j0, X.elements[i - i0], 0, j1 - j0 + 1);
}
- } catch (ArrayIndexOutOfBoundsException e) {
+ }
+ catch(ArrayIndexOutOfBoundsException e) {
throw new ArrayIndexOutOfBoundsException(ERR_REINDEX);
}
return X;
@@ -451,12 +441,13 @@ public class Matrix {
public final Matrix getMatrix(final int[] r, final int[] c) {
final Matrix X = new Matrix(r.length, c.length);
try {
- for (int i = 0; i < r.length; i++) {
- for (int j = 0; j < c.length; j++) {
+ for(int i = 0; i < r.length; i++) {
+ for(int j = 0; j < c.length; j++) {
X.elements[i][j] = elements[r[i]][c[j]];
}
}
- } catch (ArrayIndexOutOfBoundsException e) {
+ }
+ catch(ArrayIndexOutOfBoundsException e) {
throw new ArrayIndexOutOfBoundsException(ERR_REINDEX);
}
return X;
@@ -474,10 +465,11 @@ public class Matrix {
public final Matrix getMatrix(final int[] r, final int j0, final int j1) {
final Matrix X = new Matrix(r.length, j1 - j0 + 1);
try {
- for (int i = 0; i < r.length; i++) {
+ for(int i = 0; i < r.length; i++) {
System.arraycopy(elements[r[i]], j0, X.elements[i], 0, j1 - j0 + 1);
}
- } catch (ArrayIndexOutOfBoundsException e) {
+ }
+ catch(ArrayIndexOutOfBoundsException e) {
throw new ArrayIndexOutOfBoundsException(ERR_REINDEX);
}
return X;
@@ -495,12 +487,13 @@ public class Matrix {
public final Matrix getMatrix(final int i0, final int i1, final int[] c) {
final Matrix X = new Matrix(i1 - i0 + 1, c.length);
try {
- for (int i = i0; i <= i1; i++) {
- for (int j = 0; j < c.length; j++) {
+ for(int i = i0; i <= i1; i++) {
+ for(int j = 0; j < c.length; j++) {
X.elements[i - i0][j] = elements[i][c[j]];
}
}
- } catch (ArrayIndexOutOfBoundsException e) {
+ }
+ catch(ArrayIndexOutOfBoundsException e) {
throw new ArrayIndexOutOfBoundsException(ERR_REINDEX);
}
return X;
@@ -518,10 +511,11 @@ public class Matrix {
*/
public final void setMatrix(final int i0, final int i1, final int j0, final int j1, final Matrix X) {
try {
- for (int i = i0; i <= i1; i++) {
+ for(int i = i0; i <= i1; i++) {
System.arraycopy(X.elements[i - i0], 0, elements[i], j0, j1 - j0 + 1);
}
- } catch (ArrayIndexOutOfBoundsException e) {
+ }
+ catch(ArrayIndexOutOfBoundsException e) {
throw new ArrayIndexOutOfBoundsException(ERR_REINDEX);
}
}
@@ -536,12 +530,13 @@ public class Matrix {
*/
public final void setMatrix(final int[] r, final int[] c, final Matrix X) {
try {
- for (int i = 0; i < r.length; i++) {
- for (int j = 0; j < c.length; j++) {
+ for(int i = 0; i < r.length; i++) {
+ for(int j = 0; j < c.length; j++) {
elements[r[i]][c[j]] = X.elements[i][j];
}
}
- } catch (ArrayIndexOutOfBoundsException e) {
+ }
+ catch(ArrayIndexOutOfBoundsException e) {
throw new ArrayIndexOutOfBoundsException(ERR_REINDEX);
}
}
@@ -557,10 +552,11 @@ public class Matrix {
*/
public final void setMatrix(final int[] r, final int j0, final int j1, final Matrix X) {
try {
- for (int i = 0; i < r.length; i++) {
+ for(int i = 0; i < r.length; i++) {
System.arraycopy(X.elements[i], 0, elements[r[i]], j0, j1 - j0 + 1);
}
- } catch (ArrayIndexOutOfBoundsException e) {
+ }
+ catch(ArrayIndexOutOfBoundsException e) {
throw new ArrayIndexOutOfBoundsException(ERR_REINDEX);
}
}
@@ -576,12 +572,13 @@ public class Matrix {
*/
public final void setMatrix(final int i0, final int i1, final int[] c, final Matrix X) {
try {
- for (int i = i0; i <= i1; i++) {
- for (int j = 0; j < c.length; j++) {
+ for(int i = i0; i <= i1; i++) {
+ for(int j = 0; j < c.length; j++) {
elements[i][c[j]] = X.elements[i - i0][j];
}
}
- } catch (ArrayIndexOutOfBoundsException e) {
+ }
+ catch(ArrayIndexOutOfBoundsException e) {
throw new ArrayIndexOutOfBoundsException(ERR_REINDEX);
}
}
@@ -604,7 +601,7 @@ public class Matrix {
* @param row the value of the column to be set
*/
public final void setRow(final int j, final Vector row) {
- if (row.elements.length != columndimension) {
+ if(row.elements.length != columndimension) {
throw new IllegalArgumentException(ERR_MATRIX_DIMENSIONS);
}
System.arraycopy(row.elements, 0, elements[j], 0, columndimension);
@@ -618,7 +615,7 @@ public class Matrix {
*/
public final Vector getCol(final int j) {
final Vector v = new Vector(elements.length);
- for (int i = 0; i < elements.length; i++) {
+ for(int i = 0; i < elements.length; i++) {
v.elements[i] = elements[i][j];
}
return v;
@@ -631,10 +628,10 @@ public class Matrix {
* @param column the value of the column to be set
*/
public final void setCol(final int j, final Vector column) {
- if (column.elements.length != elements.length) {
+ if(column.elements.length != elements.length) {
throw new IllegalArgumentException(ERR_MATRIX_DIMENSIONS);
}
- for (int i = 0; i < elements.length; i++) {
+ for(int i = 0; i < elements.length; i++) {
elements[i][j] = column.elements[i];
}
}
@@ -646,8 +643,8 @@ public class Matrix {
*/
public final Matrix transpose() {
final Matrix X = new Matrix(columndimension, elements.length);
- for (int i = 0; i < elements.length; i++) {
- for (int j = 0; j < columndimension; j++) {
+ for(int i = 0; i < elements.length; i++) {
+ for(int j = 0; j < columndimension; j++) {
X.elements[j][i] = elements[i][j];
}
}
@@ -665,6 +662,26 @@ public class Matrix {
}
/**
+ * C = A + s
+ *
+ * @param s scalar value
+ * @return A + s in a new Matrix
+ */
+ public final Matrix plus(final double s) {
+ return copy().plusEquals(s);
+ }
+
+ /**
+ * C = A + s
+ *
+ * @param s scalar value
+ * @return A + s * E in a new Matrix
+ */
+ public final Matrix plusDiagonal(final double s) {
+ return copy().plusDiagonalEquals(s);
+ }
+
+ /**
* C = A + s * B
*
* @param B another matrix
@@ -683,8 +700,8 @@ public class Matrix {
*/
public final Matrix plusEquals(final Matrix B) {
checkMatrixDimensions(B);
- for (int i = 0; i < elements.length; i++) {
- for (int j = 0; j < columndimension; j++) {
+ for(int i = 0; i < elements.length; i++) {
+ for(int j = 0; j < columndimension; j++) {
elements[i][j] += B.elements[i][j];
}
}
@@ -692,6 +709,34 @@ public class Matrix {
}
/**
+ * A = A + s
+ *
+ * @param s constant to add to every cell
+ * @return A + s in this Matrix
+ */
+ public final Matrix plusEquals(final double s) {
+ for(int i = 0; i < elements.length; i++) {
+ for(int j = 0; j < columndimension; j++) {
+ elements[i][j] += s;
+ }
+ }
+ return this;
+ }
+
+ /**
+ * A = A + s
+ *
+ * @param s constant to add to the diagonal
+ * @return A + s in this Matrix
+ */
+ public final Matrix plusDiagonalEquals(final double s) {
+ for(int i = 0; i < elements.length && i < columndimension; i++) {
+ elements[i][i] += s;
+ }
+ return this;
+ }
+
+ /**
* A = A + s * B
*
* @param B another matrix
@@ -700,8 +745,8 @@ public class Matrix {
*/
public final Matrix plusTimesEquals(final Matrix B, final double s) {
checkMatrixDimensions(B);
- for (int i = 0; i < elements.length; i++) {
- for (int j = 0; j < columndimension; j++) {
+ for(int i = 0; i < elements.length; i++) {
+ for(int j = 0; j < columndimension; j++) {
elements[i][j] += s * B.elements[i][j];
}
}
@@ -737,8 +782,8 @@ public class Matrix {
*/
public final Matrix minusEquals(final Matrix B) {
checkMatrixDimensions(B);
- for (int i = 0; i < elements.length; i++) {
- for (int j = 0; j < columndimension; j++) {
+ for(int i = 0; i < elements.length; i++) {
+ for(int j = 0; j < columndimension; j++) {
elements[i][j] -= B.elements[i][j];
}
}
@@ -754,8 +799,8 @@ public class Matrix {
*/
public final Matrix minusTimesEquals(final Matrix B, final double s) {
checkMatrixDimensions(B);
- for (int i = 0; i < elements.length; i++) {
- for (int j = 0; j < columndimension; j++) {
+ for(int i = 0; i < elements.length; i++) {
+ for(int j = 0; j < columndimension; j++) {
elements[i][j] -= s * B.elements[i][j];
}
}
@@ -779,8 +824,8 @@ public class Matrix {
* @return replace A by s*A
*/
public final Matrix timesEquals(final double s) {
- for (int i = 0; i < elements.length; i++) {
- for (int j = 0; j < columndimension; j++) {
+ for(int i = 0; i < elements.length; i++) {
+ for(int j = 0; j < columndimension; j++) {
elements[i][j] *= s;
}
}
@@ -796,22 +841,22 @@ public class Matrix {
*/
public final Matrix times(final Matrix B) {
// Optimized implementation, exploiting the storage layout
- if (B.elements.length != this.columndimension) {
+ if(B.elements.length != this.columndimension) {
throw new IllegalArgumentException(ERR_MATRIX_INNERDIM);
}
final Matrix X = new Matrix(this.elements.length, B.columndimension);
// Optimized ala Jama. jik order.
final double[] Bcolj = new double[this.columndimension];
- for (int j = 0; j < X.columndimension; j++) {
+ for(int j = 0; j < X.columndimension; j++) {
// Make a linear copy of column j from B
- for (int k = 0; k < this.columndimension; k++) {
+ for(int k = 0; k < this.columndimension; k++) {
Bcolj[k] = B.elements[k][j];
}
// multiply it with each row from A
- for (int i = 0; i < this.elements.length; i++) {
+ for(int i = 0; i < this.elements.length; i++) {
final double[] Arowi = this.elements[i];
double s = 0;
- for (int k = 0; k < this.columndimension; k++) {
+ for(int k = 0; k < this.columndimension; k++) {
s += Arowi[k] * Bcolj[k];
}
X.elements[i][j] = s;
@@ -828,15 +873,15 @@ public class Matrix {
* @throws IllegalArgumentException Matrix inner dimensions must agree.
*/
public final Vector times(final Vector B) {
- if (B.elements.length != this.columndimension) {
+ if(B.elements.length != this.columndimension) {
throw new IllegalArgumentException(ERR_MATRIX_INNERDIM);
}
final Vector X = new Vector(this.elements.length);
// multiply it with each row from A
- for (int i = 0; i < this.elements.length; i++) {
+ for(int i = 0; i < this.elements.length; i++) {
final double[] Arowi = this.elements[i];
double s = 0;
- for (int k = 0; k < this.columndimension; k++) {
+ for(int k = 0; k < this.columndimension; k++) {
s += Arowi[k] * B.elements[k];
}
X.elements[i] = s;
@@ -852,14 +897,14 @@ public class Matrix {
* @throws IllegalArgumentException Matrix inner dimensions must agree.
*/
public final Vector transposeTimes(final Vector B) {
- if (B.elements.length != elements.length) {
+ if(B.elements.length != elements.length) {
throw new IllegalArgumentException(ERR_MATRIX_INNERDIM);
}
final Vector X = new Vector(this.columndimension);
// multiply it with each row from A
- for (int i = 0; i < this.columndimension; i++) {
+ for(int i = 0; i < this.columndimension; i++) {
double s = 0;
- for (int k = 0; k < elements.length; k++) {
+ for(int k = 0; k < elements.length; k++) {
s += elements[k][i] * B.elements[k];
}
X.elements[i] = s;
@@ -875,20 +920,20 @@ public class Matrix {
* @throws IllegalArgumentException Matrix inner dimensions must agree.
*/
public final Matrix transposeTimes(final Matrix B) {
- if (B.elements.length != elements.length) {
+ if(B.elements.length != elements.length) {
throw new IllegalArgumentException(ERR_MATRIX_INNERDIM);
}
final Matrix X = new Matrix(this.columndimension, B.columndimension);
final double[] Bcolj = new double[elements.length];
- for (int j = 0; j < X.columndimension; j++) {
+ for(int j = 0; j < X.columndimension; j++) {
// Make a linear copy of column j from B
- for (int k = 0; k < elements.length; k++) {
+ for(int k = 0; k < elements.length; k++) {
Bcolj[k] = B.elements[k][j];
}
// multiply it with each row from A
- for (int i = 0; i < this.columndimension; i++) {
+ for(int i = 0; i < this.columndimension; i++) {
double s = 0;
- for (int k = 0; k < elements.length; k++) {
+ for(int k = 0; k < elements.length; k++) {
s += elements[k][i] * Bcolj[k];
}
X.elements[i][j] = s;
@@ -905,17 +950,17 @@ public class Matrix {
* @throws IllegalArgumentException Matrix inner dimensions must agree.
*/
public final Matrix timesTranspose(final Matrix B) {
- if (B.columndimension != this.columndimension) {
+ if(B.columndimension != this.columndimension) {
throw new IllegalArgumentException(ERR_MATRIX_INNERDIM);
}
final Matrix X = new Matrix(this.elements.length, B.elements.length);
- for (int j = 0; j < X.elements.length; j++) {
+ for(int j = 0; j < X.elements.length; j++) {
final double[] Browj = B.elements[j];
// multiply it with each row from A
- for (int i = 0; i < this.elements.length; i++) {
+ for(int i = 0; i < this.elements.length; i++) {
final double[] Arowi = this.elements[i];
double s = 0;
- for (int k = 0; k < this.columndimension; k++) {
+ for(int k = 0; k < this.columndimension; k++) {
s += Arowi[k] * Browj[k];
}
X.elements[i][j] = s;
@@ -933,23 +978,23 @@ public class Matrix {
*/
public final Matrix transposeTimesTranspose(Matrix B) {
// Optimized implementation, exploiting the storage layout
- if (this.elements.length != B.columndimension) {
+ if(this.elements.length != B.columndimension) {
throw new IllegalArgumentException("Matrix inner dimensions must agree: " + getRowDimensionality() + "," + getColumnDimensionality() + " * " + B.getRowDimensionality() + "," + B.getColumnDimensionality());
}
final Matrix X = new Matrix(this.columndimension, B.elements.length);
// Optimized ala Jama. jik order.
final double[] Acolj = new double[this.elements.length];
- for (int j = 0; j < X.elements.length; j++) {
+ for(int j = 0; j < X.elements.length; j++) {
// Make a linear copy of column j from B
- for (int k = 0; k < this.elements.length; k++) {
+ for(int k = 0; k < this.elements.length; k++) {
Acolj[k] = this.elements[k][j];
}
final double[] Xrow = X.elements[j];
// multiply it with each row from A
- for (int i = 0; i < B.elements.length; i++) {
+ for(int i = 0; i < B.elements.length; i++) {
final double[] Browi = B.elements[i];
double s = 0;
- for (int k = 0; k < B.columndimension; k++) {
+ for(int k = 0; k < B.columndimension; k++) {
s += Browi[k] * Acolj[k];
}
Xrow[i] = s;
@@ -978,6 +1023,19 @@ public class Matrix {
}
/**
+ * Matrix inverse for square matrixes.
+ *
+ * @return inverse(A), or inverse(A + epsilon E) if singular.
+ */
+ public final Matrix robustInverse() {
+ LUDecomposition d = new LUDecomposition(this);
+ if(!d.isNonsingular()) {
+ d = new LUDecomposition(plusDiagonal(SINGULARITY_CHEAT).getArrayRef(), elements.length, columndimension);
+ }
+ return d.solve(identity(elements.length, elements.length));
+ }
+
+ /**
* Matrix determinant
*
* @return determinant
@@ -1011,7 +1069,7 @@ public class Matrix {
*/
public final double trace() {
double t = 0;
- for (int i = 0; i < Math.min(elements.length, columndimension); i++) {
+ for(int i = 0; i < Math.min(elements.length, columndimension); i++) {
t += elements[i][i];
}
return t;
@@ -1024,9 +1082,9 @@ public class Matrix {
*/
public final double norm1() {
double f = 0;
- for (int j = 0; j < columndimension; j++) {
+ for(int j = 0; j < columndimension; j++) {
double s = 0;
- for (int i = 0; i < elements.length; i++) {
+ for(int i = 0; i < elements.length; i++) {
s += Math.abs(elements[i][j]);
}
f = Math.max(f, s);
@@ -1050,9 +1108,9 @@ public class Matrix {
*/
public final double normInf() {
double f = 0;
- for (int i = 0; i < elements.length; i++) {
+ for(int i = 0; i < elements.length; i++) {
double s = 0;
- for (int j = 0; j < columndimension; j++) {
+ for(int j = 0; j < columndimension; j++) {
s += Math.abs(elements[i][j]);
}
f = Math.max(f, s);
@@ -1067,8 +1125,8 @@ public class Matrix {
*/
public final double normF() {
double f = 0;
- for (int i = 0; i < elements.length; i++) {
- for (int j = 0; j < columndimension; j++) {
+ for(int i = 0; i < elements.length; i++) {
+ for(int j = 0; j < columndimension; j++) {
f = MathUtil.fastHypot(f, elements[i][j]);
}
}
@@ -1079,14 +1137,14 @@ public class Matrix {
* Normalizes the columns of this matrix to length of 1.0.
*/
public final void normalizeColumns() {
- for (int col = 0; col < columndimension; col++) {
+ for(int col = 0; col < columndimension; col++) {
double norm = 0.0;
- for (int row = 0; row < elements.length; row++) {
+ for(int row = 0; row < elements.length; row++) {
norm = norm + (elements[row][col] * elements[row][col]);
}
norm = Math.sqrt(norm);
- if (norm != 0) {
- for (int row = 0; row < elements.length; row++) {
+ if(norm != 0) {
+ for(int row = 0; row < elements.length; row++) {
elements[row][col] /= norm;
}
}
@@ -1105,13 +1163,13 @@ public class Matrix {
* columns of this matrix
*/
public final boolean linearlyIndependent(final Matrix columnMatrix) {
- if (columnMatrix.columndimension != 1) {
+ if(columnMatrix.columndimension != 1) {
throw new IllegalArgumentException("a.getColumnDimension() != 1");
}
- if (this.elements.length != columnMatrix.elements.length) {
+ if(this.elements.length != columnMatrix.elements.length) {
throw new IllegalArgumentException(ERR_MATRIX_DIMENSIONS);
}
- if (this.columndimension + columnMatrix.columndimension > this.elements.length) {
+ if(this.columndimension + columnMatrix.columndimension > this.elements.length) {
return false;
}
final StringBuilder msg = LoggingConfiguration.DEBUG ? new StringBuilder() : null;
@@ -1119,20 +1177,22 @@ public class Matrix {
final double[][] a = new double[columndimension + 1][elements.length - 1];
final double[] b = new double[columndimension + 1];
- for (int i = 0; i < a.length; i++) {
- for (int j = 0; j < a[i].length; j++) {
- if (i < columndimension) {
+ for(int i = 0; i < a.length; i++) {
+ for(int j = 0; j < a[i].length; j++) {
+ if(i < columndimension) {
a[i][j] = elements[j][i];
- } else {
+ }
+ else {
a[i][j] = columnMatrix.elements[j][0];
}
}
}
- for (int i = 0; i < b.length; i++) {
- if (i < columndimension) {
+ for(int i = 0; i < b.length; i++) {
+ if(i < columndimension) {
b[i] = elements[elements.length - 1][i];
- } else {
+ }
+ else {
b[i] = columnMatrix.elements[i][0];
}
}
@@ -1143,7 +1203,7 @@ public class Matrix {
final double[][] coefficients = les.getCoefficents();
final double[] rhs = les.getRHS();
- if (msg != null) {
+ if(msg != null) {
msg.append("\na' ").append(FormatUtil.format(this.getArrayRef()));
msg.append("\nb' ").append(FormatUtil.format(columnMatrix.getColumnPackedCopy()));
@@ -1152,20 +1212,20 @@ public class Matrix {
msg.append("\nleq ").append(les.equationsToString(4));
}
- for (int i = 0; i < coefficients.length; i++) {
+ for(int i = 0; i < coefficients.length; i++) {
boolean allCoefficientsZero = true;
- for (int j = 0; j < coefficients[i].length; j++) {
+ for(int j = 0; j < coefficients[i].length; j++) {
final double value = coefficients[i][j];
- if (Math.abs(value) > DELTA) {
+ if(Math.abs(value) > DELTA) {
allCoefficientsZero = false;
break;
}
}
// allCoefficients=0 && rhs=0 -> linearly dependent
- if (allCoefficientsZero) {
+ if(allCoefficientsZero) {
final double value = rhs[i];
- if (Math.abs(value) < DELTA) {
- if (msg != null) {
+ if(Math.abs(value) < DELTA) {
+ if(msg != null) {
msg.append("\nvalue ").append(value).append('[').append(i).append(']');
msg.append("\nlinearly independent ").append(false);
Logger.getLogger(this.getClass().getName()).fine(msg.toString());
@@ -1175,7 +1235,7 @@ public class Matrix {
}
}
- if (msg != null) {
+ if(msg != null) {
msg.append("\nlinearly independent ").append(true);
Logger.getLogger(this.getClass().getName()).fine(msg.toString());
}
@@ -1183,128 +1243,17 @@ public class Matrix {
}
/**
- * Returns a matrix derived by Gauss-Jordan-elimination using RationalNumbers
- * for the transformations.
- *
- * @return a matrix derived by Gauss-Jordan-elimination using RationalNumbers
- * for the transformations
- */
- public final Matrix exactGaussJordanElimination() {
- final RationalNumber[][] gauss = exactGaussElimination();
-
- // reduced form
- for (int row = gauss.length - 1; row > 0; row--) {
- int firstCol = -1;
- for (int col = 0; col < gauss[row].length && firstCol == -1; col++) {
- // if(gauss.get(row, col) != 0.0) // i.e. == 1
- if (gauss[row][col].equals(RationalNumber.ONE)) {
- firstCol = col;
- }
- }
- if (firstCol > -1) {
- for (int currentRow = row - 1; currentRow >= 0; currentRow--) {
- RationalNumber multiplier = gauss[currentRow][firstCol].copy();
- for (int col = firstCol; col < gauss[currentRow].length; col++) {
- RationalNumber subtrahent = gauss[row][col].times(multiplier);
- gauss[currentRow][col] = gauss[currentRow][col].minus(subtrahent);
- }
- }
- }
- }
- return new Matrix(gauss);
- }
-
- /**
- * Perform an exact Gauss-elimination of this Matrix using RationalNumbers to
- * yield highest possible accuracy.
- *
- * @return an array of arrays of RationalNumbers representing the
- * Gauss-eliminated form of this Matrix
- */
- private final RationalNumber[][] exactGaussElimination() {
- final RationalNumber[][] gauss = new RationalNumber[elements.length][this.columndimension];
- for (int row = 0; row < elements.length; row++) {
- for (int col = 0; col < this.columndimension; col++) {
- gauss[row][col] = new RationalNumber(elements[row][col]);
- }
- }
- return exactGaussElimination(gauss);
- }
-
- /**
- * Perform recursive Gauss-elimination on the given matrix of RationalNumbers.
- *
- * @param gauss an array of arrays of RationalNumber
- * @return recursive derived Gauss-elimination-form of the given matrix of
- * RationalNumbers
- */
- private static final RationalNumber[][] exactGaussElimination(final RationalNumber[][] gauss) {
- int firstCol = -1;
- int firstRow = -1;
-
- // 1. find first column unequal to zero
- for (int col = 0; col < gauss[0].length && firstCol == -1; col++) {
- for (int row = 0; row < gauss.length && firstCol == -1; row++) {
- // if(gauss.get(row, col) != 0.0)
- if (!gauss[row][col].equals(RationalNumber.ZERO)) {
- firstCol = col;
- firstRow = row;
- }
- }
- }
-
- // 2. set row as first row
- if (firstCol != -1) {
- if (firstRow != 0) {
- final RationalNumber[] row = new RationalNumber[gauss[firstRow].length];
- System.arraycopy(gauss[firstRow], 0, row, 0, gauss[firstRow].length);
- System.arraycopy(gauss[0], 0, gauss[firstRow], 0, gauss[firstRow].length);
- System.arraycopy(row, 0, gauss[0], 0, row.length);
- }
-
- // 3. create leading 1
- if (!gauss[0][firstCol].equals(RationalNumber.ONE)) {
- final RationalNumber inverse = gauss[0][firstCol].multiplicativeInverse();
- for (int col = 0; col < gauss[0].length; col++) {
- gauss[0][col] = gauss[0][col].times(inverse);
- }
- }
-
- // 4. eliminate values unequal to zero below leading 1
- for (int row = 1; row < gauss.length; row++) {
- final RationalNumber multiplier = gauss[row][firstCol].copy();
- // if(multiplier != 0.0)
- if (!multiplier.equals(RationalNumber.ZERO)) {
- for (int col = firstCol; col < gauss[row].length; col++) {
- final RationalNumber subtrahent = gauss[0][col].times(multiplier);
- gauss[row][col] = gauss[row][col].minus(subtrahent);
- }
- }
- }
-
- // 5. recursion
- if (gauss.length > 1) {
- final RationalNumber[][] subMatrix = new RationalNumber[gauss.length - 1][gauss[1].length];
- System.arraycopy(gauss, 1, subMatrix, 0, gauss.length - 1);
- final RationalNumber[][] eliminatedSubMatrix = exactGaussElimination(subMatrix);
- System.arraycopy(eliminatedSubMatrix, 0, gauss, 1, eliminatedSubMatrix.length);
- }
- }
- return gauss;
- }
-
- /**
* Returns true, if this matrix is symmetric, false otherwise.
*
* @return true, if this matrix is symmetric, false otherwise
*/
public final boolean isSymmetric() {
- if (elements.length != columndimension) {
+ if(elements.length != columndimension) {
return false;
}
- for (int i = 0; i < elements.length; i++) {
- for (int j = i + 1; j < columndimension; j++) {
- if (elements[i][j] != elements[j][i]) {
+ for(int i = 0; i < elements.length; i++) {
+ for(int j = i + 1; j < columndimension; j++) {
+ if(elements[i][j] != elements[j][i]) {
return false;
}
}
@@ -1321,16 +1270,17 @@ public class Matrix {
public final Matrix completeBasis() {
Matrix basis = copy();
Matrix result = null;
- for (int i = 0; i < elements.length; i++) {
+ for(int i = 0; i < elements.length; i++) {
final Matrix e_i = new Matrix(elements.length, 1);
e_i.elements[0][i] = 1.0;
final boolean li = basis.linearlyIndependent(e_i);
// TODO: efficiency - appendColumns is expensive.
- if (li) {
- if (result == null) {
+ if(li) {
+ if(result == null) {
result = e_i.copy();
- } else {
+ }
+ else {
result = result.appendColumns(e_i);
}
basis = basis.appendColumns(e_i);
@@ -1348,16 +1298,17 @@ public class Matrix {
public final Matrix completeToOrthonormalBasis() {
Matrix basis = copy();
Matrix result = null;
- for (int i = 0; i < elements.length; i++) {
+ for(int i = 0; i < elements.length; i++) {
final Matrix e_i = new Matrix(elements.length, 1);
e_i.elements[i][0] = 1.0;
final boolean li = basis.linearlyIndependent(e_i);
// TODO: efficiency - appendColumns is expensive.
- if (li) {
- if (result == null) {
+ if(li) {
+ if(result == null) {
result = e_i.copy();
- } else {
+ }
+ else {
result = result.appendColumns(e_i);
}
basis = basis.appendColumns(e_i);
@@ -1374,16 +1325,17 @@ public class Matrix {
* @return the new matrix with the appended columns
*/
public final Matrix appendColumns(final Matrix columns) {
- if (elements.length != columns.elements.length) {
+ if(elements.length != columns.elements.length) {
throw new IllegalArgumentException(ERR_MATRIX_DIMENSIONS);
}
final Matrix result = new Matrix(elements.length, columndimension + columns.columndimension);
- for (int i = 0; i < result.columndimension; i++) {
+ for(int i = 0; i < result.columndimension; i++) {
// FIXME: optimize - excess copying!
- if (i < columndimension) {
+ if(i < columndimension) {
result.setCol(i, getCol(i));
- } else {
+ }
+ else {
result.setCol(i, columns.getCol(i - columndimension));
}
}
@@ -1399,10 +1351,10 @@ public class Matrix {
Matrix v = copy();
// FIXME: optimize - excess copying!
- for (int i = 1; i < columndimension; i++) {
+ for(int i = 1; i < columndimension; i++) {
final Vector u_i = getCol(i);
final Vector sum = new Vector(elements.length);
- for (int j = 0; j < i; j++) {
+ for(int j = 0; j < i; j++) {
final Vector v_j = v.getCol(j);
double scalar = u_i.transposeTimes(v_j) / v_j.transposeTimes(v_j);
sum.plusTimesEquals(v_j, scalar);
@@ -1425,7 +1377,7 @@ public class Matrix {
*/
public final Matrix cheatToAvoidSingularity(final double constant) {
final Matrix a = this.copy();
- for (int i = 0; i < a.columndimension && i < a.elements.length; i++) {
+ for(int i = 0; i < a.columndimension && i < a.elements.length; i++) {
// if(a.get(i, i) < constant)
{
a.elements[i][i] += constant;
@@ -1460,40 +1412,40 @@ public class Matrix {
TDoubleArrayList v = new TDoubleArrayList();
// Ignore initial empty lines
- while (tokenizer.nextToken() == StreamTokenizer.TT_EOL) {
+ while(tokenizer.nextToken() == StreamTokenizer.TT_EOL) {
// ignore initial empty lines
}
- if (tokenizer.ttype == StreamTokenizer.TT_EOF) {
+ if(tokenizer.ttype == StreamTokenizer.TT_EOF) {
throw new java.io.IOException("Unexpected EOF on matrix read.");
}
do {
v.add(FormatUtil.parseDouble(tokenizer.sval)); // Read & store 1st
// row.
}
- while (tokenizer.nextToken() == StreamTokenizer.TT_WORD);
+ while(tokenizer.nextToken() == StreamTokenizer.TT_WORD);
int n = v.size(); // Now we've got the number of columns!
double row[] = v.toArray();
ArrayList<double[]> rowV = new ArrayList<>();
rowV.add(row); // Start storing rows instead of columns.
- while (tokenizer.nextToken() == StreamTokenizer.TT_WORD) {
+ while(tokenizer.nextToken() == StreamTokenizer.TT_WORD) {
// While non-empty lines
rowV.add(row = new double[n]);
int j = 0;
do {
- if (j >= n) {
+ if(j >= n) {
throw new java.io.IOException("Row " + v.size() + " is too long.");
}
row[j++] = FormatUtil.parseDouble(tokenizer.sval);
}
- while (tokenizer.nextToken() == StreamTokenizer.TT_WORD);
- if (j < n) {
+ while(tokenizer.nextToken() == StreamTokenizer.TT_WORD);
+ if(j < n) {
throw new java.io.IOException("Row " + v.size() + " is too short.");
}
}
int m = rowV.size(); // Now we've got the number of rows.
double[][] A = new double[m][];
- for (int i = 0; i < m; i++) {
+ for(int i = 0; i < m; i++) {
A[i] = rowV.get(i);
}
return new Matrix(A);
@@ -1503,7 +1455,7 @@ public class Matrix {
* Check if size(A) == size(B)
*/
protected void checkMatrixDimensions(Matrix B) {
- if (B.getRowDimensionality() != getRowDimensionality() || B.getColumnDimensionality() != getColumnDimensionality()) {
+ if(B.getRowDimensionality() != getRowDimensionality() || B.getColumnDimensionality() != getColumnDimensionality()) {
throw new IllegalArgumentException("Matrix dimensions must agree.");
}
}
@@ -1520,25 +1472,25 @@ public class Matrix {
@Override
public boolean equals(Object obj) {
- if (this == obj) {
+ if(this == obj) {
return true;
}
- if (obj == null) {
+ if(obj == null) {
return false;
}
- if (getClass() != obj.getClass()) {
+ if(getClass() != obj.getClass()) {
return false;
}
final Matrix other = (Matrix) obj;
- if (this.elements.length != other.elements.length) {
+ if(this.elements.length != other.elements.length) {
return false;
}
- if (this.columndimension != other.columndimension) {
+ if(this.columndimension != other.columndimension) {
return false;
}
- for (int i = 0; i < this.elements.length; i++) {
- for (int j = 0; j < this.columndimension; j++) {
- if (this.elements[i][j] != other.elements[i][j]) {
+ for(int i = 0; i < this.elements.length; i++) {
+ for(int j = 0; j < this.columndimension; j++) {
+ if(this.elements[i][j] != other.elements[i][j]) {
return false;
}
}
@@ -1555,25 +1507,25 @@ public class Matrix {
* @return true if delta smaller than maximum
*/
public boolean almostEquals(Object obj, double maxdelta) {
- if (this == obj) {
+ if(this == obj) {
return true;
}
- if (obj == null) {
+ if(obj == null) {
return false;
}
- if (getClass() != obj.getClass()) {
+ if(getClass() != obj.getClass()) {
return false;
}
final Matrix other = (Matrix) obj;
- if (this.elements.length != other.elements.length) {
+ if(this.elements.length != other.elements.length) {
return false;
}
- if (this.columndimension != other.columndimension) {
+ if(this.columndimension != other.columndimension) {
return false;
}
- for (int i = 0; i < this.elements.length; i++) {
- for (int j = 0; j < this.columndimension; j++) {
- if (Math.abs(this.elements[i][j] - other.elements[i][j]) > maxdelta) {
+ for(int i = 0; i < this.elements.length; i++) {
+ for(int j = 0; j < this.columndimension; j++) {
+ if(Math.abs(this.elements[i][j] - other.elements[i][j]) > maxdelta) {
return false;
}
}
generated by cgit v1.2.3 (git 2.25.1) at 2024-05-21 21:30:06 +0000