0001: /* AUTO-GENERATED */
0002: package JSci.maths.matrices;
0003:
0004: import JSci.GlobalSettings;
0005: import JSci.maths.ArrayMath;
0006: import JSci.maths.ExtraMath;
0007: import JSci.maths.LinearMath;
0008: import JSci.maths.MaximumIterationsExceededException;
0009: import JSci.maths.Mapping;
0010: import JSci.maths.vectors.AbstractDoubleVector;
0011: import JSci.maths.groups.AbelianGroup;
0012: import JSci.maths.algebras.*;
0013: import JSci.maths.fields.*;
0014:
0015: /**
0016: * The AbstractDoubleSquareMatrix class provides an object for encapsulating double square matrices.
0017: * @version 2.2
0018: * @author Mark Hale
0019: */
0020: public abstract class AbstractDoubleSquareMatrix extends
0021: AbstractDoubleMatrix implements SquareMatrix {
0022: protected transient AbstractDoubleSquareMatrix[] LU;
0023: protected transient int[] LUpivot;
0024:
0025: /**
0026: * Constructs a matrix.
0027: */
0028: protected AbstractDoubleSquareMatrix(final int size) {
0029: super (size, size);
0030: }
0031:
0032: /**
0033: * Converts this matrix to an integer matrix.
0034: * @return an integer square matrix
0035: */
0036: public AbstractIntegerMatrix toIntegerMatrix() {
0037: final int ans[][] = new int[numRows][numCols];
0038: for (int i = 0; i < numRows; i++) {
0039: for (int j = 0; j < numCols; j++)
0040: ans[i][j] = Math.round((float) getElement(i, j));
0041: }
0042: return new IntegerSquareMatrix(ans);
0043: }
0044:
0045: /**
0046: * Converts this matrix to a complex matrix.
0047: * @return a complex square matrix
0048: */
0049: public AbstractComplexMatrix toComplexMatrix() {
0050: ComplexSquareMatrix cm = new ComplexSquareMatrix(numRows);
0051: for (int i = 0; i < numRows; i++) {
0052: for (int j = 0; j < numCols; j++)
0053: cm.setElement(i, j, getElement(i, j), 0.0);
0054: }
0055: return cm;
0056: }
0057:
0058: /**
0059: * Returns true if this matrix is symmetric.
0060: */
0061: public boolean isSymmetric() {
0062: return this .equals(this .transpose());
0063: }
0064:
0065: /**
0066: * Returns true if this matrix is unitary.
0067: */
0068: public boolean isUnitary() {
0069: return this .multiply(this .transpose()).equals(
0070: DoubleDiagonalMatrix.identity(numRows));
0071: }
0072:
0073: /**
0074: * Returns the determinant.
0075: */
0076: public double det() {
0077: if (numRows == 2) {
0078: return getElement(0, 0) * getElement(1, 1)
0079: - getElement(0, 1) * getElement(1, 0);
0080: } else {
0081: final AbstractDoubleSquareMatrix lu[] = this
0082: .luDecompose(null);
0083: double det = lu[1].getElement(0, 0);
0084: for (int i = 1; i < numRows; i++)
0085: det *= lu[1].getElement(i, i);
0086: return det * LUpivot[numRows];
0087: }
0088: }
0089:
0090: /**
0091: * Returns the trace.
0092: */
0093: public double trace() {
0094: double result = getElement(0, 0);
0095: for (int i = 1; i < numRows; i++)
0096: result += getElement(i, i);
0097: return result;
0098: }
0099:
0100: /**
0101: * Returns the operator norm.
0102: * @exception MaximumIterationsExceededException If it takes more than 50 iterations to determine an eigenvalue.
0103: */
0104: public double operatorNorm()
0105: throws MaximumIterationsExceededException {
0106: return Math
0107: .sqrt(ArrayMath
0108: .max(LinearMath
0109: .eigenvalueSolveSymmetric((AbstractDoubleSquareMatrix) (this
0110: .transpose().multiply(this )))));
0111: }
0112:
0113: //============
0114: // OPERATIONS
0115: //============
0116:
0117: /**
0118: * Returns the negative of this matrix.
0119: */
0120: public AbelianGroup.Member negate() {
0121: final double array[][] = new double[numRows][numCols];
0122: for (int i = 0; i < numRows; i++) {
0123: array[i][0] = -getElement(i, 0);
0124: for (int j = 1; j < numCols; j++)
0125: array[i][j] = -getElement(i, j);
0126: }
0127: return new DoubleSquareMatrix(array);
0128: }
0129:
0130: // ADDITION
0131:
0132: /**
0133: * Returns the addition of this matrix and another.
0134: * @param m a double square matrix
0135: * @exception MatrixDimensionException If the matrices are not square or different sizes.
0136: */
0137: public final AbstractDoubleMatrix add(final AbstractDoubleMatrix m) {
0138: if (m instanceof AbstractDoubleSquareMatrix)
0139: return add((AbstractDoubleSquareMatrix) m);
0140: else if (numRows == m.rows() && numCols == m.columns())
0141: return add(new SquareMatrixAdaptor(m));
0142: else
0143: throw new MatrixDimensionException(
0144: "Matrices are different sizes.");
0145: }
0146:
0147: /**
0148: * Returns the addition of this matrix and another.
0149: * @param m a double square matrix
0150: * @exception MatrixDimensionException If the matrices are different sizes.
0151: */
0152: public AbstractDoubleSquareMatrix add(
0153: final AbstractDoubleSquareMatrix m) {
0154: if (numRows == m.rows() && numCols == m.columns()) {
0155: final double array[][] = new double[numRows][numCols];
0156: for (int i = 0; i < numRows; i++) {
0157: array[i][0] = getElement(i, 0) + m.getElement(i, 0);
0158: for (int j = 1; j < numCols; j++)
0159: array[i][j] = getElement(i, j) + m.getElement(i, j);
0160: }
0161: return new DoubleSquareMatrix(array);
0162: } else {
0163: throw new MatrixDimensionException(
0164: "Matrices are different sizes.");
0165: }
0166: }
0167:
0168: // SUBTRACTION
0169:
0170: /**
0171: * Returns the subtraction of this matrix and another.
0172: * @param m a double square matrix
0173: * @exception MatrixDimensionException If the matrices are not square or different sizes.
0174: */
0175: public final AbstractDoubleMatrix subtract(
0176: final AbstractDoubleMatrix m) {
0177: if (m instanceof AbstractDoubleSquareMatrix)
0178: return subtract((AbstractDoubleSquareMatrix) m);
0179: else if (numRows == m.rows() && numCols == m.columns())
0180: return subtract(new SquareMatrixAdaptor(m));
0181: else
0182: throw new MatrixDimensionException(
0183: "Matrices are different sizes.");
0184: }
0185:
0186: /**
0187: * Returns the subtraction of this matrix by another.
0188: * @param m a double square matrix
0189: * @exception MatrixDimensionException If the matrices are different sizes.
0190: */
0191: public AbstractDoubleSquareMatrix subtract(
0192: final AbstractDoubleSquareMatrix m) {
0193: if (numRows == m.rows() && numCols == m.columns()) {
0194: final double array[][] = new double[numRows][numCols];
0195: for (int i = 0; i < numRows; i++) {
0196: array[i][0] = getElement(i, 0) - m.getElement(i, 0);
0197: for (int j = 1; j < numCols; j++)
0198: array[i][j] = getElement(i, j) - m.getElement(i, j);
0199: }
0200: return new DoubleSquareMatrix(array);
0201: } else {
0202: throw new MatrixDimensionException(
0203: "Matrices are different sizes.");
0204: }
0205: }
0206:
0207: // SCALAR MULTIPLICATION
0208:
0209: /**
0210: * Returns the multiplication of this matrix by a scalar.
0211: * @param x a double.
0212: * @return a double square matrix.
0213: */
0214: public AbstractDoubleMatrix scalarMultiply(final double x) {
0215: final double array[][] = new double[numRows][numCols];
0216: for (int i = 0; i < numRows; i++) {
0217: array[i][0] = x * getElement(i, 0);
0218: for (int j = 1; j < numCols; j++)
0219: array[i][j] = x * getElement(i, j);
0220: }
0221: return new DoubleSquareMatrix(array);
0222: }
0223:
0224: // SCALAR DIVISON
0225:
0226: /**
0227: * Returns the division of this matrix by a scalar.
0228: * @param x a double.
0229: * @return a double square matrix.
0230: */
0231: public AbstractDoubleMatrix scalarDivide(final double x) {
0232: final double array[][] = new double[numRows][numCols];
0233: for (int i = 0; i < numRows; i++) {
0234: array[i][0] = getElement(i, 0) / x;
0235: for (int j = 1; j < numCols; j++)
0236: array[i][j] = getElement(i, j) / x;
0237: }
0238: return new DoubleSquareMatrix(array);
0239: }
0240:
0241: // SCALAR PRODUCT
0242:
0243: /**
0244: * Returns the scalar product of this matrix and another.
0245: * @param m a double square matrix.
0246: * @exception MatrixDimensionException If the matrices are not square or different sizes.
0247: */
0248: public final double scalarProduct(final AbstractDoubleMatrix m) {
0249: if (m instanceof AbstractDoubleSquareMatrix)
0250: return scalarProduct((AbstractDoubleSquareMatrix) m);
0251: else if (numRows == m.rows() && numCols == m.columns())
0252: return scalarProduct(new SquareMatrixAdaptor(m));
0253: else
0254: throw new MatrixDimensionException(
0255: "Matrices are different sizes.");
0256: }
0257:
0258: /**
0259: * Returns the scalar product of this matrix and another.
0260: * @param m a double square matrix.
0261: * @exception MatrixDimensionException If the matrices are different sizes.
0262: */
0263: public double scalarProduct(final AbstractDoubleSquareMatrix m) {
0264: if (numRows == m.rows() && numCols == m.columns()) {
0265: double ans = 0;
0266: for (int i = 0; i < numRows; i++) {
0267: ans += getElement(i, 0) * m.getElement(i, 0);
0268: for (int j = 1; j < numCols; j++)
0269: ans += getElement(i, j) * m.getElement(i, j);
0270: }
0271: return ans;
0272: } else {
0273: throw new MatrixDimensionException(
0274: "Matrices are different sizes.");
0275: }
0276: }
0277:
0278: // MATRIX MULTIPLICATION
0279:
0280: /**
0281: * Returns the multiplication of this matrix and another.
0282: * @param m a double square matrix
0283: * @exception MatrixDimensionException If the matrices are different sizes.
0284: */
0285: public AbstractDoubleSquareMatrix multiply(
0286: final AbstractDoubleSquareMatrix m) {
0287: if (numCols == m.rows()) {
0288: final int mColumns = m.columns();
0289: final double array[][] = new double[numRows][mColumns];
0290: for (int j = 0; j < numRows; j++) {
0291: for (int k = 0; k < mColumns; k++) {
0292: array[j][k] = getElement(j, 0) * m.getElement(0, k);
0293: for (int n = 1; n < numCols; n++)
0294: array[j][k] += getElement(j, n)
0295: * m.getElement(n, k);
0296: }
0297: }
0298: return new DoubleSquareMatrix(array);
0299: } else {
0300: throw new MatrixDimensionException("Incompatible matrices.");
0301: }
0302: }
0303:
0304: // DIRECT SUM
0305:
0306: /**
0307: * Returns the direct sum of this matrix and another.
0308: */
0309: public AbstractDoubleSquareMatrix directSum(
0310: final AbstractDoubleSquareMatrix m) {
0311: final double array[][] = new double[numRows + m.numRows][numCols
0312: + m.numCols];
0313: for (int i = 0; i < numRows; i++) {
0314: for (int j = 0; j < numCols; j++)
0315: array[i][j] = getElement(i, j);
0316: }
0317: for (int i = 0; i < m.numRows; i++) {
0318: for (int j = 0; j < m.numCols; j++)
0319: array[i + numRows][j + numCols] = m.getElement(i, j);
0320: }
0321: return new DoubleSquareMatrix(array);
0322: }
0323:
0324: // TENSOR PRODUCT
0325:
0326: /**
0327: * Returns the tensor product of this matrix and another.
0328: */
0329: public AbstractDoubleSquareMatrix tensor(
0330: final AbstractDoubleSquareMatrix m) {
0331: final double array[][] = new double[numRows * m.numRows][numCols
0332: * m.numCols];
0333: for (int i = 0; i < numRows; i++) {
0334: for (int j = 0; j < numCols; j++) {
0335: for (int k = 0; k < m.numRows; j++) {
0336: for (int l = 0; l < m.numCols; l++)
0337: array[i * m.numRows + k][j * m.numCols + l] = getElement(
0338: i, j)
0339: * m.getElement(k, l);
0340: }
0341: }
0342: }
0343: return new DoubleSquareMatrix(array);
0344: }
0345:
0346: // TRANSPOSE
0347:
0348: /**
0349: * Returns the transpose of this matrix.
0350: * @return a double square matrix
0351: */
0352: public Matrix transpose() {
0353: final double array[][] = new double[numCols][numRows];
0354: for (int i = 0; i < numRows; i++) {
0355: array[0][i] = getElement(i, 0);
0356: for (int j = 1; j < numCols; j++)
0357: array[j][i] = getElement(i, j);
0358: }
0359: return new DoubleSquareMatrix(array);
0360: }
0361:
0362: // INVERSE
0363:
0364: /**
0365: * Returns the inverse of this matrix.
0366: * @return a double square matrix.
0367: */
0368: public AbstractDoubleSquareMatrix inverse() {
0369: final int N = numRows;
0370: final double arrayL[][] = new double[N][N];
0371: final double arrayU[][] = new double[N][N];
0372: final AbstractDoubleSquareMatrix lu[] = this .luDecompose(null);
0373: arrayL[0][0] = 1.0 / lu[0].getElement(0, 0);
0374: arrayU[0][0] = 1.0 / lu[1].getElement(0, 0);
0375: for (int i = 1; i < N; i++) {
0376: arrayL[i][i] = 1.0 / lu[0].getElement(i, i);
0377: arrayU[i][i] = 1.0 / lu[1].getElement(i, i);
0378: }
0379: for (int i = 0; i < N - 1; i++) {
0380: for (int j = i + 1; j < N; j++) {
0381: double tmpL = 0.0, tmpU = 0.0;
0382: for (int k = i; k < j; k++) {
0383: tmpL -= lu[0].getElement(j, k) * arrayL[k][i];
0384: tmpU -= arrayU[i][k] * lu[1].getElement(k, j);
0385: }
0386: arrayL[j][i] = tmpL / lu[0].getElement(j, j);
0387: arrayU[i][j] = tmpU / lu[1].getElement(j, j);
0388: }
0389: }
0390: // matrix multiply arrayU x arrayL
0391: final double inv[][] = new double[N][N];
0392: for (int i = 0; i < N; i++) {
0393: for (int j = 0; j < i; j++) {
0394: for (int k = i; k < N; k++)
0395: inv[i][LUpivot[j]] += arrayU[i][k] * arrayL[k][j];
0396: }
0397: for (int j = i; j < N; j++) {
0398: for (int k = j; k < N; k++)
0399: inv[i][LUpivot[j]] += arrayU[i][k] * arrayL[k][j];
0400: }
0401: }
0402: return new DoubleSquareMatrix(inv);
0403: }
0404:
0405: // LU DECOMPOSITION
0406:
0407: /**
0408: * Returns the LU decomposition of this matrix.
0409: * @param pivot an empty array of length <code>rows()+1</code>
0410: * to hold the pivot information (null if not interested).
0411: * The last array element will contain the parity.
0412: * @return an array with [0] containing the L-matrix
0413: * and [1] containing the U-matrix.
0414: * @jsci.planetmath LUDecomposition
0415: */
0416: public AbstractDoubleSquareMatrix[] luDecompose(int pivot[]) {
0417: if (LU != null) {
0418: if (pivot != null)
0419: System.arraycopy(LUpivot, 0, pivot, 0, pivot.length);
0420: return LU;
0421: }
0422: int pivotrow;
0423: final int N = numRows;
0424: final double arrayL[][] = new double[N][N];
0425: final double arrayU[][] = new double[N][N];
0426: if (pivot == null)
0427: pivot = new int[N + 1];
0428: for (int i = 0; i < N; i++)
0429: pivot[i] = i;
0430: pivot[N] = 1;
0431: // LU decomposition to arrayU
0432: for (int j = 0; j < N; j++) {
0433: for (int i = 0; i < j; i++) {
0434: double tmp = getElement(pivot[i], j);
0435: for (int k = 0; k < i; k++)
0436: tmp -= arrayU[i][k] * arrayU[k][j];
0437: arrayU[i][j] = tmp;
0438: }
0439: double max = 0.0;
0440: pivotrow = j;
0441: for (int i = j; i < N; i++) {
0442: double tmp = getElement(pivot[i], j);
0443: for (int k = 0; k < j; k++)
0444: tmp -= arrayU[i][k] * arrayU[k][j];
0445: arrayU[i][j] = tmp;
0446: // while we're here search for a pivot for arrayU[j][j]
0447: tmp = Math.abs(tmp);
0448: if (tmp > max) {
0449: max = tmp;
0450: pivotrow = i;
0451: }
0452: }
0453: // swap row j with pivotrow
0454: if (pivotrow != j) {
0455: double[] tmprow = arrayU[j];
0456: arrayU[j] = arrayU[pivotrow];
0457: arrayU[pivotrow] = tmprow;
0458: int k = pivot[j];
0459: pivot[j] = pivot[pivotrow];
0460: pivot[pivotrow] = k;
0461: // update parity
0462: pivot[N] = -pivot[N];
0463: }
0464: // divide by pivot
0465: double tmp = arrayU[j][j];
0466: for (int i = j + 1; i < N; i++)
0467: arrayU[i][j] /= tmp;
0468: }
0469: // move lower triangular part to arrayL
0470: for (int j = 0; j < N; j++) {
0471: arrayL[j][j] = 1.0;
0472: for (int i = j + 1; i < N; i++) {
0473: arrayL[i][j] = arrayU[i][j];
0474: arrayU[i][j] = 0.0;
0475: }
0476: }
0477: LU = new AbstractDoubleSquareMatrix[2];
0478: LU[0] = new DoubleSquareMatrix(arrayL);
0479: LU[1] = new DoubleSquareMatrix(arrayU);
0480: LUpivot = new int[pivot.length];
0481: System.arraycopy(pivot, 0, LUpivot, 0, pivot.length);
0482: return LU;
0483: }
0484:
0485: /**
0486: * Returns the LU decomposition of this matrix.
0487: * Warning: no pivoting.
0488: * @return an array with [0] containing the L-matrix
0489: * and [1] containing the U-matrix.
0490: * @jsci.planetmath LUDecomposition
0491: */
0492: public AbstractDoubleSquareMatrix[] luDecompose() {
0493: final int N = numRows;
0494: final double arrayL[][] = new double[N][N];
0495: final double arrayU[][] = new double[N][N];
0496: // LU decomposition to arrayU
0497: for (int j = 0; j < N; j++) {
0498: for (int i = 0; i < j; i++) {
0499: double tmp = getElement(i, j);
0500: for (int k = 0; k < i; k++)
0501: tmp -= arrayU[i][k] * arrayU[k][j];
0502: arrayU[i][j] = tmp;
0503: }
0504: for (int i = j; i < N; i++) {
0505: double tmp = getElement(i, j);
0506: for (int k = 0; k < j; k++)
0507: tmp -= arrayU[i][k] * arrayU[k][j];
0508: arrayU[i][j] = tmp;
0509: }
0510: // divide
0511: double tmp = arrayU[j][j];
0512: for (int i = j + 1; i < N; i++)
0513: arrayU[i][j] /= tmp;
0514: }
0515: // move lower triangular part to arrayL
0516: for (int j = 0; j < N; j++) {
0517: arrayL[j][j] = 1.0;
0518: for (int i = j + 1; i < N; i++) {
0519: arrayL[i][j] = arrayU[i][j];
0520: arrayU[i][j] = 0.0;
0521: }
0522: }
0523: AbstractDoubleSquareMatrix[] lu = new AbstractDoubleSquareMatrix[2];
0524: lu[0] = new DoubleSquareMatrix(arrayL);
0525: lu[1] = new DoubleSquareMatrix(arrayU);
0526: return lu;
0527: }
0528:
0529: // CHOLESKY DECOMPOSITION
0530:
0531: /**
0532: * Returns the Cholesky decomposition of this matrix.
0533: * Matrix must be symmetric and positive definite.
0534: * @return an array with [0] containing the L-matrix and [1] containing the U-matrix.
0535: * @jsci.planetmath CholeskyDecomposition
0536: */
0537: public AbstractDoubleSquareMatrix[] choleskyDecompose() {
0538: final int N = numRows;
0539: final double arrayL[][] = new double[N][N];
0540: final double arrayU[][] = new double[N][N];
0541: double tmp = Math.sqrt(getElement(0, 0));
0542: arrayL[0][0] = arrayU[0][0] = tmp;
0543: for (int i = 1; i < N; i++)
0544: arrayL[i][0] = arrayU[0][i] = getElement(i, 0) / tmp;
0545: for (int j = 1; j < N; j++) {
0546: tmp = getElement(j, j);
0547: for (int i = 0; i < j; i++)
0548: tmp -= arrayL[j][i] * arrayL[j][i];
0549: arrayL[j][j] = arrayU[j][j] = Math.sqrt(tmp);
0550: for (int i = j + 1; i < N; i++) {
0551: tmp = getElement(i, j);
0552: for (int k = 0; k < i; k++)
0553: tmp -= arrayL[j][k] * arrayU[k][i];
0554: arrayL[i][j] = arrayU[j][i] = tmp / arrayU[j][j];
0555: }
0556: }
0557: final AbstractDoubleSquareMatrix lu[] = new AbstractDoubleSquareMatrix[2];
0558: lu[0] = new DoubleSquareMatrix(arrayL);
0559: lu[1] = new DoubleSquareMatrix(arrayU);
0560: return lu;
0561: }
0562:
0563: // QR DECOMPOSITION
0564:
0565: /**
0566: * Returns the QR decomposition of this matrix.
0567: * Based on the code from <a href="http://math.nist.gov/javanumerics/jama/">JAMA</a> (public domain).
0568: * @return an array with [0] containing the Q-matrix and [1] containing the R-matrix.
0569: * @jsci.planetmath QRDecomposition
0570: */
0571: public AbstractDoubleSquareMatrix[] qrDecompose() {
0572: final int N = numRows;
0573: final double array[][] = new double[N][N];
0574: final double arrayQ[][] = new double[N][N];
0575: final double arrayR[][] = new double[N][N];
0576: // copy matrix
0577: for (int i = 0; i < N; i++) {
0578: array[i][0] = getElement(i, 0);
0579: for (int j = 1; j < N; j++)
0580: array[i][j] = getElement(i, j);
0581: }
0582: for (int k = 0; k < N; k++) {
0583: // compute l2-norm of kth column
0584: double norm = array[k][k];
0585: for (int i = k + 1; i < N; i++)
0586: norm = ExtraMath.hypot(norm, array[i][k]);
0587: if (norm != 0.0) {
0588: // form kth Householder vector
0589: if (array[k][k] < 0.0)
0590: norm = -norm;
0591: for (int i = k; i < N; i++)
0592: array[i][k] /= norm;
0593: array[k][k] += 1.0;
0594: // apply transformation to remaining columns
0595: for (int j = k + 1; j < N; j++) {
0596: double s = array[k][k] * array[k][j];
0597: for (int i = k + 1; i < N; i++)
0598: s += array[i][k] * array[i][j];
0599: s /= array[k][k];
0600: for (int i = k; i < N; i++)
0601: array[i][j] -= s * array[i][k];
0602: }
0603: }
0604: arrayR[k][k] = -norm;
0605: }
0606: for (int k = N - 1; k >= 0; k--) {
0607: arrayQ[k][k] = 1.0;
0608: for (int j = k; j < N; j++) {
0609: if (array[k][k] != 0.0) {
0610: double s = array[k][k] * arrayQ[k][j];
0611: for (int i = k + 1; i < N; i++)
0612: s += array[i][k] * arrayQ[i][j];
0613: s /= array[k][k];
0614: for (int i = k; i < N; i++)
0615: arrayQ[i][j] -= s * array[i][k];
0616: }
0617: }
0618: }
0619: for (int i = 0; i < N; i++) {
0620: for (int j = i + 1; j < N; j++)
0621: arrayR[i][j] = array[i][j];
0622: }
0623: final AbstractDoubleSquareMatrix qr[] = new AbstractDoubleSquareMatrix[2];
0624: qr[0] = new DoubleSquareMatrix(arrayQ);
0625: qr[1] = new DoubleSquareMatrix(arrayR);
0626: return qr;
0627: }
0628:
0629: // SINGULAR VALUE DECOMPOSITION
0630:
0631: /**
0632: * Returns the singular value decomposition of this matrix.
0633: * Based on the code from <a href="http://math.nist.gov/javanumerics/jama/">JAMA</a> (public domain).
0634: * @return an array with [0] containing the U-matrix, [1] containing the S-matrix and [2] containing the V-matrix.
0635: * @jsci.planetmath SingularValueDecomposition
0636: */
0637: public AbstractDoubleSquareMatrix[] singularValueDecompose() {
0638: final int N = numRows;
0639: final int Nm1 = N - 1;
0640: final double array[][] = new double[N][N];
0641: final double arrayU[][] = new double[N][N];
0642: final double arrayS[] = new double[N];
0643: final double arrayV[][] = new double[N][N];
0644: final double e[] = new double[N];
0645: final double work[] = new double[N];
0646: // copy matrix
0647: for (int i = 0; i < N; i++) {
0648: array[i][0] = getElement(i, 0);
0649: for (int j = 1; j < N; j++)
0650: array[i][j] = getElement(i, j);
0651: }
0652: // reduce matrix to bidiagonal form
0653: for (int k = 0; k < Nm1; k++) {
0654: // compute the transformation for the kth column
0655: // compute l2-norm of kth column
0656: arrayS[k] = array[k][k];
0657: for (int i = k + 1; i < N; i++)
0658: arrayS[k] = ExtraMath.hypot(arrayS[k], array[i][k]);
0659: if (arrayS[k] != 0.0) {
0660: if (array[k][k] < 0.0)
0661: arrayS[k] = -arrayS[k];
0662: for (int i = k; i < N; i++)
0663: array[i][k] /= arrayS[k];
0664: array[k][k] += 1.0;
0665: }
0666: arrayS[k] = -arrayS[k];
0667: for (int j = k + 1; j < N; j++) {
0668: if (arrayS[k] != 0.0) {
0669: // apply the transformation
0670: double t = array[k][k] * array[k][j];
0671: for (int i = k + 1; i < N; i++)
0672: t += array[i][k] * array[i][j];
0673: t /= array[k][k];
0674: for (int i = k; i < N; i++)
0675: array[i][j] -= t * array[i][k];
0676: }
0677: e[j] = array[k][j];
0678: }
0679: for (int i = k; i < N; i++)
0680: arrayU[i][k] = array[i][k];
0681: if (k < N - 2) {
0682: // compute the kth row transformation
0683: // compute l2-norm of kth column
0684: e[k] = e[k + 1];
0685: for (int i = k + 2; i < N; i++)
0686: e[k] = ExtraMath.hypot(e[k], e[i]);
0687: if (e[k] != 0.0) {
0688: if (e[k + 1] < 0.0)
0689: e[k] = -e[k];
0690: for (int i = k + 1; i < N; i++)
0691: e[i] /= e[k];
0692: e[k + 1] += 1.0;
0693: }
0694: e[k] = -e[k];
0695: if (e[k] != 0.0) {
0696: // apply the transformation
0697: for (int i = k + 1; i < N; i++) {
0698: work[i] = 0.0;
0699: for (int j = k + 1; j < N; j++)
0700: work[i] += e[j] * array[i][j];
0701: }
0702: for (int j = k + 1; j < N; j++) {
0703: double t = e[j] / e[k + 1];
0704: for (int i = k + 1; i < N; i++)
0705: array[i][j] -= t * work[i];
0706: }
0707: }
0708: for (int i = k + 1; i < N; i++)
0709: arrayV[i][k] = e[i];
0710: }
0711: }
0712: // setup the final bidiagonal matrix of order p
0713: int p = N;
0714: arrayS[Nm1] = array[Nm1][Nm1];
0715: e[N - 2] = array[N - 2][Nm1];
0716: e[Nm1] = 0.0;
0717: arrayU[Nm1][Nm1] = 1.0;
0718: for (int k = N - 2; k >= 0; k--) {
0719: if (arrayS[k] != 0.0) {
0720: for (int j = k + 1; j < N; j++) {
0721: double t = arrayU[k][k] * arrayU[k][j];
0722: for (int i = k + 1; i < N; i++)
0723: t += arrayU[i][k] * arrayU[i][j];
0724: t /= arrayU[k][k];
0725: for (int i = k; i < N; i++)
0726: arrayU[i][j] -= t * arrayU[i][k];
0727: }
0728: for (int i = k; i < N; i++)
0729: arrayU[i][k] = -arrayU[i][k];
0730: arrayU[k][k] += 1.0;
0731: for (int i = 0; i < k - 1; i++)
0732: arrayU[i][k] = 0.0;
0733: } else {
0734: for (int i = 0; i < N; i++)
0735: arrayU[i][k] = 0.0;
0736: arrayU[k][k] = 1.0;
0737: }
0738: }
0739: for (int k = Nm1; k >= 0; k--) {
0740: if (k < N - 2 && e[k] != 0.0) {
0741: for (int j = k + 1; j < N; j++) {
0742: double t = arrayV[k + 1][k] * arrayV[k + 1][j];
0743: for (int i = k + 2; i < N; i++)
0744: t += arrayV[i][k] * arrayV[i][j];
0745: t /= arrayV[k + 1][k];
0746: for (int i = k + 1; i < N; i++)
0747: arrayV[i][j] -= t * arrayV[i][k];
0748: }
0749: }
0750: for (int i = 0; i < N; i++)
0751: arrayV[i][k] = 0.0;
0752: arrayV[k][k] = 1.0;
0753: }
0754: final double eps = Math.pow(2.0, -52.0);
0755: int iter = 0;
0756: while (p > 0) {
0757: int k, action;
0758: // action = 1 if arrayS[p] and e[k-1] are negligible and k<p
0759: // action = 2 if arrayS[k] is negligible and k<p
0760: // action = 3 if e[k-1] is negligible, k<p, and arrayS[k], ..., arrayS[p] are not negligible (QR step)
0761: // action = 4 if e[p-1] is negligible (convergence)
0762: for (k = p - 2; k >= -1; k--) {
0763: if (k == -1)
0764: break;
0765: if (Math.abs(e[k]) <= eps
0766: * (Math.abs(arrayS[k]) + Math
0767: .abs(arrayS[k + 1]))) {
0768: e[k] = 0.0;
0769: break;
0770: }
0771: }
0772: if (k == p - 2) {
0773: action = 4;
0774: } else {
0775: int ks;
0776: for (ks = p - 1; ks >= k; ks--) {
0777: if (ks == k)
0778: break;
0779: double t = (ks != p ? Math.abs(e[ks]) : 0.0)
0780: + (ks != k + 1 ? Math.abs(e[ks - 1]) : 0.0);
0781: if (Math.abs(arrayS[ks]) <= eps * t) {
0782: arrayS[ks] = 0.0;
0783: break;
0784: }
0785: }
0786: if (ks == k) {
0787: action = 3;
0788: } else if (ks == p - 1) {
0789: action = 1;
0790: } else {
0791: action = 2;
0792: k = ks;
0793: }
0794: }
0795: k++;
0796: switch (action) {
0797: // deflate negligible arrayS[p]
0798: case 1: {
0799: double f = e[p - 2];
0800: e[p - 2] = 0.0;
0801: for (int j = p - 2; j >= k; j--) {
0802: double t = ExtraMath.hypot(arrayS[j], f);
0803: final double cs = arrayS[j] / t;
0804: final double sn = f / t;
0805: arrayS[j] = t;
0806: if (j != k) {
0807: f = -sn * e[j - 1];
0808: e[j - 1] *= cs;
0809: }
0810: for (int i = 0; i < N; i++) {
0811: t = cs * arrayV[i][j] + sn * arrayV[i][p - 1];
0812: arrayV[i][p - 1] = -sn * arrayV[i][j] + cs
0813: * arrayV[i][p - 1];
0814: arrayV[i][j] = t;
0815: }
0816: }
0817: }
0818: break;
0819: // split at negligible arrayS[k]
0820: case 2: {
0821: double f = e[k - 1];
0822: e[k - 1] = 0.0;
0823: for (int j = k; j < p; j++) {
0824: double t = ExtraMath.hypot(arrayS[j], f);
0825: final double cs = arrayS[j] / t;
0826: final double sn = f / t;
0827: arrayS[j] = t;
0828: f = -sn * e[j];
0829: e[j] *= cs;
0830: for (int i = 0; i < N; i++) {
0831: t = cs * arrayU[i][j] + sn * arrayU[i][k - 1];
0832: arrayU[i][k - 1] = -sn * arrayU[i][j] + cs
0833: * arrayU[i][k - 1];
0834: arrayU[i][j] = t;
0835: }
0836: }
0837: }
0838: break;
0839: // perform one QR step
0840: case 3: {
0841: // calculate the shift
0842: final double scale = Math.max(Math.max(Math.max(Math
0843: .max(Math.abs(arrayS[p - 1]), Math
0844: .abs(arrayS[p - 2])), Math
0845: .abs(e[p - 2])), Math.abs(arrayS[k])), Math
0846: .abs(e[k]));
0847: double sp = arrayS[p - 1] / scale;
0848: double spm1 = arrayS[p - 2] / scale;
0849: double epm1 = e[p - 2] / scale;
0850: double sk = arrayS[k] / scale;
0851: double ek = e[k] / scale;
0852: double b = ((spm1 + sp) * (spm1 - sp) + epm1 * epm1) / 2.0;
0853: double c = (sp * epm1) * (sp * epm1);
0854: double shift = 0.0;
0855: if (b != 0.0 || c != 0.0) {
0856: shift = Math.sqrt(b * b + c);
0857: if (b < 0.0)
0858: shift = -shift;
0859: shift = c / (b + shift);
0860: }
0861: double f = (sk + sp) * (sk - sp) + shift;
0862: double g = sk * ek;
0863: // chase zeros
0864: for (int j = k; j < p - 1; j++) {
0865: double t = ExtraMath.hypot(f, g);
0866: double cs = f / t;
0867: double sn = g / t;
0868: if (j != k)
0869: e[j - 1] = t;
0870: f = cs * arrayS[j] + sn * e[j];
0871: e[j] = cs * e[j] - sn * arrayS[j];
0872: g = sn * arrayS[j + 1];
0873: arrayS[j + 1] *= cs;
0874: for (int i = 0; i < N; i++) {
0875: t = cs * arrayV[i][j] + sn * arrayV[i][j + 1];
0876: arrayV[i][j + 1] = -sn * arrayV[i][j] + cs
0877: * arrayV[i][j + 1];
0878: arrayV[i][j] = t;
0879: }
0880: t = ExtraMath.hypot(f, g);
0881: cs = f / t;
0882: sn = g / t;
0883: arrayS[j] = t;
0884: f = cs * e[j] + sn * arrayS[j + 1];
0885: arrayS[j + 1] = -sn * e[j] + cs * arrayS[j + 1];
0886: g = sn * e[j + 1];
0887: e[j + 1] *= cs;
0888: if (j < Nm1) {
0889: for (int i = 0; i < N; i++) {
0890: t = cs * arrayU[i][j] + sn
0891: * arrayU[i][j + 1];
0892: arrayU[i][j + 1] = -sn * arrayU[i][j] + cs
0893: * arrayU[i][j + 1];
0894: arrayU[i][j] = t;
0895: }
0896: }
0897: }
0898: e[p - 2] = f;
0899: iter++;
0900: }
0901: break;
0902: // convergence
0903: case 4: {
0904: // make the singular values positive
0905: if (arrayS[k] <= 0.0) {
0906: arrayS[k] = -arrayS[k];
0907: for (int i = 0; i < p; i++)
0908: arrayV[i][k] = -arrayV[i][k];
0909: }
0910: // order the singular values
0911: while (k < p - 1) {
0912: if (arrayS[k] >= arrayS[k + 1])
0913: break;
0914: double tmp = arrayS[k];
0915: arrayS[k] = arrayS[k + 1];
0916: arrayS[k + 1] = tmp;
0917: if (k < Nm1) {
0918: for (int i = 0; i < N; i++) {
0919: tmp = arrayU[i][k + 1];
0920: arrayU[i][k + 1] = arrayU[i][k];
0921: arrayU[i][k] = tmp;
0922: tmp = arrayV[i][k + 1];
0923: arrayV[i][k + 1] = arrayV[i][k];
0924: arrayV[i][k] = tmp;
0925: }
0926: }
0927: k++;
0928: }
0929: iter = 0;
0930: p--;
0931: }
0932: break;
0933: }
0934: }
0935: final AbstractDoubleSquareMatrix svd[] = new AbstractDoubleSquareMatrix[3];
0936: svd[0] = new DoubleSquareMatrix(arrayU);
0937: svd[1] = new DoubleDiagonalMatrix(arrayS);
0938: svd[2] = new DoubleSquareMatrix(arrayV);
0939: return svd;
0940: }
0941:
0942: // POLAR DECOMPOSITION
0943:
0944: /**
0945: * Returns the polar decomposition of this matrix.
0946: */
0947: public AbstractDoubleSquareMatrix[] polarDecompose() {
0948: final int N = numRows;
0949: final AbstractDoubleVector evec[] = new AbstractDoubleVector[N];
0950: double eval[];
0951: try {
0952: eval = LinearMath.eigenSolveSymmetric(this , evec);
0953: } catch (MaximumIterationsExceededException e) {
0954: return null;
0955: }
0956: final double tmpa[][] = new double[N][N];
0957: final double tmpm[][] = new double[N][N];
0958: double abs;
0959: for (int i = 0; i < N; i++) {
0960: abs = Math.abs(eval[i]);
0961: tmpa[i][0] = eval[i] * evec[i].getComponent(0) / abs;
0962: tmpm[i][0] = abs * evec[i].getComponent(0);
0963: for (int j = 1; j < N; j++) {
0964: tmpa[i][j] = eval[i] * evec[i].getComponent(j) / abs;
0965: tmpm[i][j] = abs * evec[i].getComponent(j);
0966: }
0967: }
0968: final double arg[][] = new double[N][N];
0969: final double mod[][] = new double[N][N];
0970: for (int i = 0; i < N; i++) {
0971: for (int j = 0; j < N; j++) {
0972: arg[i][j] = evec[0].getComponent(i) * tmpa[0][j];
0973: mod[i][j] = evec[0].getComponent(i) * tmpm[0][j];
0974: for (int k = 1; k < N; k++) {
0975: arg[i][j] += evec[k].getComponent(i) * tmpa[k][j];
0976: mod[i][j] += evec[k].getComponent(i) * tmpm[k][j];
0977: }
0978: }
0979: }
0980: final AbstractDoubleSquareMatrix us[] = new AbstractDoubleSquareMatrix[2];
0981: us[0] = new DoubleSquareMatrix(arg);
0982: us[1] = new DoubleSquareMatrix(mod);
0983: return us;
0984: }
0985:
0986: // MAP ELEMENTS
0987:
0988: /**
0989: * Applies a function on all the matrix elements.
0990: * @param f a user-defined function
0991: * @return a double square matrix
0992: */
0993: public AbstractDoubleMatrix mapElements(final Mapping f) {
0994: final double array[][] = new double[numRows][numCols];
0995: for (int i = 0; i < numRows; i++) {
0996: array[i][0] = f.map(getElement(i, 0));
0997: for (int j = 1; j < numCols; j++)
0998: array[i][j] = f.map(getElement(i, j));
0999: }
1000: return new DoubleSquareMatrix(array);
1001: }
1002:
1003: private static class SquareMatrixAdaptor extends
1004: AbstractDoubleSquareMatrix {
1005: private final AbstractDoubleMatrix matrix;
1006:
1007: private SquareMatrixAdaptor(AbstractDoubleMatrix m) {
1008: super (m.rows());
1009: matrix = m;
1010: }
1011:
1012: public double getElement(int i, int j) {
1013: return matrix.getElement(i, j);
1014: }
1015:
1016: public void setElement(int i, int j, double x) {
1017: matrix.setElement(i, j, x);
1018: }
1019: }
1020: }
|