Class Algebra

java.lang.Object
cern.colt.PersistentObject
cern.colt.matrix.linalg.Algebra
All Implemented Interfaces:
Serializable, Cloneable

public class Algebra extends PersistentObject
Linear algebraic matrix operations operating on DoubleMatrix2D; concentrates most functionality of this package.
Version:
1.0, 09/24/99
See Also:
  • Field Details

    • DEFAULT

      public static final Algebra DEFAULT
      A default Algebra object; has Property.DEFAULT attached for tolerance. Allows ommiting to construct an Algebra object time and again. Note that this Algebra object is immutable. Any attempt to assign a new Property object to it (via method setProperty), or to alter the tolerance of its property object (via property().setTolerance(...)) will throw an exception.
    • ZERO

      public static final Algebra ZERO
      A default Algebra object; has Property.ZERO attached for tolerance. Allows ommiting to construct an Algebra object time and again. Note that this Algebra object is immutable. Any attempt to assign a new Property object to it (via method setProperty), or to alter the tolerance of its property object (via property().setTolerance(...)) will throw an exception.
    • property

      protected Property property
      The property object attached to this instance.
  • Constructor Details

    • Algebra

      public Algebra()
      Constructs a new instance with an equality tolerance given by Property.DEFAULT.tolerance().
    • Algebra

      public Algebra(double tolerance)
      Constructs a new instance with the given equality tolerance.
      Parameters:
      tolerance - the tolerance to be used for equality operations.
  • Method Details

    • chol

      private CholeskyDecomposition chol(DoubleMatrix2D matrix)
      Constructs and returns the cholesky-decomposition of the given matrix.
    • clone

      public Object clone()
      Returns a copy of the receiver. The attached property object is also copied. Hence, the property object of the copy is mutable.
      Overrides:
      clone in class PersistentObject
      Returns:
      a copy of the receiver.
    • cond

      public double cond(DoubleMatrix2D A)
      Returns the condition of matrix A, which is the ratio of largest to smallest singular value.
    • det

      public double det(DoubleMatrix2D A)
      Returns the determinant of matrix A.
      Returns:
      the determinant.
    • eig

      Constructs and returns the Eigenvalue-decomposition of the given matrix.
    • hypot

      protected static double hypot(double a, double b)
      Returns sqrt(a^2 + b^2) without under/overflow.
    • hypotFunction

      protected static DoubleDoubleFunction hypotFunction()
      Returns sqrt(a^2 + b^2) without under/overflow.
    • inverse

      public DoubleMatrix2D inverse(DoubleMatrix2D A)
      Returns the inverse or pseudo-inverse of matrix A.
      Returns:
      a new independent matrix; inverse(matrix) if the matrix is square, pseudoinverse otherwise.
    • lu

      private LUDecomposition lu(DoubleMatrix2D matrix)
      Constructs and returns the LU-decomposition of the given matrix.
    • mult

      public double mult(DoubleMatrix1D x, DoubleMatrix1D y)
      Inner product of two vectors; Sum(x[i] * y[i]). Also known as dot product.
      Equivalent to x.zDotProduct(y).
      Parameters:
      x - the first source vector.
      y - the second source matrix.
      Returns:
      the inner product.
      Throws:
      IllegalArgumentException - if x.size() != y.size().
    • mult

      Linear algebraic matrix-vector multiplication; z = A * y. z[i] = Sum(A[i,j] * y[j]), i=0..A.rows()-1, j=0..y.size()-1.
      Parameters:
      A - the source matrix.
      y - the source vector.
      Returns:
      z; a new vector with z.size()==A.rows().
      Throws:
      IllegalArgumentException - if A.columns() != y.size().
    • mult

      Linear algebraic matrix-matrix multiplication; C = A x B. C[i,j] = Sum(A[i,k] * B[k,j]), k=0..n-1.
      Matrix shapes: A(m x n), B(n x p), C(m x p).
      Parameters:
      A - the first source matrix.
      B - the second source matrix.
      Returns:
      C; a new matrix holding the results, with C.rows()=A.rows(), C.columns()==B.columns().
      Throws:
      IllegalArgumentException - if B.rows() != A.columns().
    • multOuter

      Outer product of two vectors; Sets A[i,j] = x[i] * y[j].
      Parameters:
      x - the first source vector.
      y - the second source vector.
      A - the matrix to hold the results. Set this parameter to null to indicate that a new result matrix shall be constructed.
      Returns:
      A (for convenience only).
      Throws:
      IllegalArgumentException - if A.rows() != x.size() || A.columns() != y.size().
    • norm1

      public double norm1(DoubleMatrix1D x)
      Returns the one-norm of vector x, which is Sum(abs(x[i])).
    • norm1

      public double norm1(DoubleMatrix2D A)
      Returns the one-norm of matrix A, which is the maximum absolute column sum.
    • norm2

      public double norm2(DoubleMatrix1D x)
      Returns the two-norm (aka euclidean norm) of vector x; equivalent to mult(x,x).
    • norm2

      public double norm2(DoubleMatrix2D A)
      Returns the two-norm of matrix A, which is the maximum singular value; obtained from SVD.
    • normF

      public double normF(DoubleMatrix2D A)
      Returns the Frobenius norm of matrix A, which is Sqrt(Sum(A[i,j]2)).
    • normInfinity

      public double normInfinity(DoubleMatrix1D x)
      Returns the infinity norm of vector x, which is Max(abs(x[i])).
    • normInfinity

      public double normInfinity(DoubleMatrix2D A)
      Returns the infinity norm of matrix A, which is the maximum absolute row sum.
    • permute

      public DoubleMatrix1D permute(DoubleMatrix1D A, int[] indexes, double[] work)
      Modifies the given vector A such that it is permuted as specified; Useful for pivoting. Cell A[i] will go into cell A[indexes[i]].

      Example:

      Reordering
      [A,B,C,D,E] with indexes [0,4,2,3,1] yields 
      [A,E,C,D,B]
      In other words A[0]invalid input: '<'--A[0], A[1]invalid input: '<'--A[4], A[2]invalid input: '<'--A[2], A[3]invalid input: '<'--A[3], A[4]invalid input: '<'--A[1].
      
      Reordering
      [A,B,C,D,E] with indexes [0,4,1,2,3] yields 
      [A,E,B,C,D]
      In other words A[0]invalid input: '<'--A[0], A[1]invalid input: '<'--A[4], A[2]invalid input: '<'--A[1], A[3]invalid input: '<'--A[2], A[4]invalid input: '<'--A[3].
      
      Parameters:
      A - the vector to permute.
      indexes - the permutation indexes, must satisfy indexes.length==A.size() invalid input: '&'invalid input: '&' indexes[i] >= 0 invalid input: '&'invalid input: '&' indexes[i] invalid input: '<' A.size();
      work - the working storage, must satisfy work.length >= A.size(); set work==null if you don't care about performance.
      Returns:
      the modified A (for convenience only).
      Throws:
      IndexOutOfBoundsException - if indexes.length != A.size().
    • permute

      public DoubleMatrix2D permute(DoubleMatrix2D A, int[] rowIndexes, int[] columnIndexes)
      Constructs and returns a new row and column permuted selection view of matrix A; equivalent to DoubleMatrix2D.viewSelection(int[],int[]). The returned matrix is backed by this matrix, so changes in the returned matrix are reflected in this matrix, and vice-versa. Use idioms like result = permute(...).copy() to generate an independent sub matrix.
      Returns:
      the new permuted selection view.
    • permuteColumns

      public DoubleMatrix2D permuteColumns(DoubleMatrix2D A, int[] indexes, int[] work)
      Modifies the given matrix A such that it's columns are permuted as specified; Useful for pivoting. Column A[i] will go into column A[indexes[i]]. Equivalent to permuteRows(transpose(A), indexes, work).
      Parameters:
      A - the matrix to permute.
      indexes - the permutation indexes, must satisfy indexes.length==A.columns() invalid input: '&'invalid input: '&' indexes[i] >= 0 invalid input: '&'invalid input: '&' indexes[i] invalid input: '<' A.columns();
      work - the working storage, must satisfy work.length >= A.columns(); set work==null if you don't care about performance.
      Returns:
      the modified A (for convenience only).
      Throws:
      IndexOutOfBoundsException - if indexes.length != A.columns().
    • permuteRows

      public DoubleMatrix2D permuteRows(DoubleMatrix2D A, int[] indexes, int[] work)
      Modifies the given matrix A such that it's rows are permuted as specified; Useful for pivoting. Row A[i] will go into row A[indexes[i]].

      Example:

      Reordering
      [A,B,C,D,E] with indexes [0,4,2,3,1] yields 
      [A,E,C,D,B]
      In other words A[0]invalid input: '<'--A[0], A[1]invalid input: '<'--A[4], A[2]invalid input: '<'--A[2], A[3]invalid input: '<'--A[3], A[4]invalid input: '<'--A[1].
      
      Reordering
      [A,B,C,D,E] with indexes [0,4,1,2,3] yields 
      [A,E,B,C,D]
      In other words A[0]invalid input: '<'--A[0], A[1]invalid input: '<'--A[4], A[2]invalid input: '<'--A[1], A[3]invalid input: '<'--A[2], A[4]invalid input: '<'--A[3].
      
      Parameters:
      A - the matrix to permute.
      indexes - the permutation indexes, must satisfy indexes.length==A.rows() invalid input: '&'invalid input: '&' indexes[i] >= 0 invalid input: '&'invalid input: '&' indexes[i] invalid input: '<' A.rows();
      work - the working storage, must satisfy work.length >= A.rows(); set work==null if you don't care about performance.
      Returns:
      the modified A (for convenience only).
      Throws:
      IndexOutOfBoundsException - if indexes.length != A.rows().
    • pow

      public DoubleMatrix2D pow(DoubleMatrix2D A, int p)
      Linear algebraic matrix power; B = Ak invalid input: '<'==> B = A*A*...*A.
      • p >= 1: B = A*A*...*A.
      • p == 0: B = identity matrix.
      • p < 0: B = pow(inverse(A),-p).
      Implementation: Based on logarithms of 2, memory usage minimized.
      Parameters:
      A - the source matrix; must be square; stays unaffected by this operation.
      p - the exponent, can be any number.
      Returns:
      B, a newly constructed result matrix; storage-independent of A.
      Throws:
      IllegalArgumentException - if !property().isSquare(A).
    • property

      public Property property()
      Returns the property object attached to this Algebra, defining tolerance.
      Returns:
      the Property object.
      See Also:
    • qr

      private QRDecomposition qr(DoubleMatrix2D matrix)
      Constructs and returns the QR-decomposition of the given matrix.
    • rank

      public int rank(DoubleMatrix2D A)
      Returns the effective numerical rank of matrix A, obtained from Singular Value Decomposition.
    • setProperty

      public void setProperty(Property property)
      Attaches the given property object to this Algebra, defining tolerance.
      Parameters:
      the - Property object to be attached.
      Throws:
      UnsupportedOperationException - if this==DEFAULT invalid input: '&'invalid input: '&' property!=this.property() - The DEFAULT Algebra object is immutable.
      UnsupportedOperationException - if this==ZERO invalid input: '&'invalid input: '&' property!=this.property() - The ZERO Algebra object is immutable.
      See Also:
    • solve

      Solves A*X = B.
      Returns:
      X; a new independent matrix; solution if A is square, least squares solution otherwise.
    • solveTranspose

      public DoubleMatrix2D solveTranspose(DoubleMatrix2D A, DoubleMatrix2D B)
      Solves X*A = B, which is also A'*X' = B'.
      Returns:
      X; a new independent matrix; solution if A is square, least squares solution otherwise.
    • subMatrix

      private DoubleMatrix2D subMatrix(DoubleMatrix2D A, int[] rowIndexes, int columnFrom, int columnTo)
      Copies the columns of the indicated rows into a new sub matrix. sub[0..rowIndexes.length-1,0..columnTo-columnFrom] = A[rowIndexes(:),columnFrom..columnTo]; The returned matrix is not backed by this matrix, so changes in the returned matrix are not reflected in this matrix, and vice-versa.
      Parameters:
      A - the source matrix to copy from.
      rowIndexes - the indexes of the rows to copy. May be unsorted.
      columnFrom - the index of the first column to copy (inclusive).
      columnTo - the index of the last column to copy (inclusive).
      Returns:
      a new sub matrix; with sub.rows()==rowIndexes.length; sub.columns()==columnTo-columnFrom+1.
      Throws:
      IndexOutOfBoundsException - if columnFrominvalid input: '<'0 || columnTo-columnFrom+1invalid input: '<'0 || columnTo+1>matrix.columns() || for any row=rowIndexes[i]: row invalid input: '<' 0 || row >= matrix.rows().
    • subMatrix

      private DoubleMatrix2D subMatrix(DoubleMatrix2D A, int rowFrom, int rowTo, int[] columnIndexes)
      Copies the rows of the indicated columns into a new sub matrix. sub[0..rowTo-rowFrom,0..columnIndexes.length-1] = A[rowFrom..rowTo,columnIndexes(:)]; The returned matrix is not backed by this matrix, so changes in the returned matrix are not reflected in this matrix, and vice-versa.
      Parameters:
      A - the source matrix to copy from.
      rowFrom - the index of the first row to copy (inclusive).
      rowTo - the index of the last row to copy (inclusive).
      columnIndexes - the indexes of the columns to copy. May be unsorted.
      Returns:
      a new sub matrix; with sub.rows()==rowTo-rowFrom+1; sub.columns()==columnIndexes.length.
      Throws:
      IndexOutOfBoundsException - if rowFrominvalid input: '<'0 || rowTo-rowFrom+1invalid input: '<'0 || rowTo+1>matrix.rows() || for any col=columnIndexes[i]: col invalid input: '<' 0 || col >= matrix.columns().
    • subMatrix

      public DoubleMatrix2D subMatrix(DoubleMatrix2D A, int fromRow, int toRow, int fromColumn, int toColumn)
      Constructs and returns a new sub-range view which is the sub matrix A[fromRow..toRow,fromColumn..toColumn]. The returned matrix is backed by this matrix, so changes in the returned matrix are reflected in this matrix, and vice-versa. Use idioms like result = subMatrix(...).copy() to generate an independent sub matrix.
      Parameters:
      A - the source matrix.
      fromRow - The index of the first row (inclusive).
      toRow - The index of the last row (inclusive).
      fromColumn - The index of the first column (inclusive).
      toColumn - The index of the last column (inclusive).
      Returns:
      a new sub-range view.
      Throws:
      IndexOutOfBoundsException - if fromColumninvalid input: '<'0 || toColumn-fromColumn+1invalid input: '<'0 || toColumn>=A.columns() || fromRowinvalid input: '<'0 || toRow-fromRow+1invalid input: '<'0 || toRow>=A.rows()
    • svd

      Constructs and returns the SingularValue-decomposition of the given matrix.
    • toString

      public String toString(DoubleMatrix2D matrix)
      Returns a String with (propertyName, propertyValue) pairs. Useful for debugging or to quickly get the rough picture. For example,
      cond          : 14.073264490042144
      det           : Illegal operation or error: Matrix must be square.
      norm1         : 0.9620244354009628
      norm2         : 3.0
      normF         : 1.304841791648992
      normInfinity  : 1.5406551198102534
      rank          : 3
      trace         : 0
      
    • toVerboseString

      public String toVerboseString(DoubleMatrix2D matrix)
      Returns the results of toString(A) and additionally the results of all sorts of decompositions applied to the given matrix. Useful for debugging or to quickly get the rough picture. For example,
      A = 3 x 3 matrix
      249  66  68
      104 214 108
      144 146 293
      
      cond         : 3.931600417472078
      det          : 9638870.0
      norm1        : 497.0
      norm2        : 473.34508217011404
      normF        : 516.873292016525
      normInfinity : 583.0
      rank         : 3
      trace        : 756.0
      
      density                      : 1.0
      isDiagonal                   : false
      isDiagonallyDominantByColumn : true
      isDiagonallyDominantByRow    : true
      isIdentity                   : false
      isLowerBidiagonal            : false
      isLowerTriangular            : false
      isNonNegative                : true
      isOrthogonal                 : false
      isPositive                   : true
      isSingular                   : false
      isSkewSymmetric              : false
      isSquare                     : true
      isStrictlyLowerTriangular    : false
      isStrictlyTriangular         : false
      isStrictlyUpperTriangular    : false
      isSymmetric                  : false
      isTriangular                 : false
      isTridiagonal                : false
      isUnitTriangular             : false
      isUpperBidiagonal            : false
      isUpperTriangular            : false
      isZero                       : false
      lowerBandwidth               : 2
      semiBandwidth                : 3
      upperBandwidth               : 2
      
      -----------------------------------------------------------------------------
      LUDecompositionQuick(A) --> isNonSingular(A), det(A), pivot, L, U, inverse(A)
      -----------------------------------------------------------------------------
      isNonSingular = true
      det = 9638870.0
      pivot = [0, 1, 2]
      
      L = 3 x 3 matrix
      1        0       0
      0.417671 1       0
      0.578313 0.57839 1
      
      U = 3 x 3 matrix
      249  66         68       
        0 186.433735  79.598394
        0   0        207.635819
      
      inverse(A) = 3 x 3 matrix
       0.004869 -0.000976 -0.00077 
      -0.001548  0.006553 -0.002056
      -0.001622 -0.002786  0.004816
      
      -----------------------------------------------------------------
      QRDecomposition(A) --> hasFullRank(A), H, Q, R, pseudo inverse(A)
      -----------------------------------------------------------------
      hasFullRank = true
      
      H = 3 x 3 matrix
      1.814086 0        0
      0.34002  1.903675 0
      0.470797 0.428218 2
      
      Q = 3 x 3 matrix
      -0.814086  0.508871  0.279845
      -0.34002  -0.808296  0.48067 
      -0.470797 -0.296154 -0.831049
      
      R = 3 x 3 matrix
      -305.864349 -195.230337 -230.023539
         0        -182.628353  467.703164
         0           0        -309.13388 
      
      pseudo inverse(A) = 3 x 3 matrix
       0.006601  0.001998 -0.005912
      -0.005105  0.000444  0.008506
      -0.000905 -0.001555  0.002688
      
      --------------------------------------------------------------------------
      CholeskyDecomposition(A) --> isSymmetricPositiveDefinite(A), L, inverse(A)
      --------------------------------------------------------------------------
      isSymmetricPositiveDefinite = false
      
      L = 3 x 3 matrix
      15.779734  0         0       
       6.590732 13.059948  0       
       9.125629  6.573948 12.903724
      
      inverse(A) = Illegal operation or error: Matrix is not symmetric positive definite.
      
      ---------------------------------------------------------------------
      EigenvalueDecomposition(A) --> D, V, realEigenvalues, imagEigenvalues
      ---------------------------------------------------------------------
      realEigenvalues = 1 x 3 matrix
      462.796507 172.382058 120.821435
      imagEigenvalues = 1 x 3 matrix
      0 0 0
      
      D = 3 x 3 matrix
      462.796507   0          0       
        0        172.382058   0       
        0          0        120.821435
      
      V = 3 x 3 matrix
      -0.398877 -0.778282  0.094294
      -0.500327  0.217793 -0.806319
      -0.768485  0.66553   0.604862
      
      ---------------------------------------------------------------------
      SingularValueDecomposition(A) --> cond(A), rank(A), norm2(A), U, S, V
      ---------------------------------------------------------------------
      cond = 3.931600417472078
      rank = 3
      norm2 = 473.34508217011404
      
      U = 3 x 3 matrix
      0.46657  -0.877519  0.110777
      0.50486   0.161382 -0.847982
      0.726243  0.45157   0.51832 
      
      S = 3 x 3 matrix
      473.345082   0          0       
        0        169.137441   0       
        0          0        120.395013
      
      V = 3 x 3 matrix
      0.577296 -0.808174  0.116546
      0.517308  0.251562 -0.817991
      0.631761  0.532513  0.563301
      
    • trace

      public double trace(DoubleMatrix2D A)
      Returns the sum of the diagonal elements of matrix A; Sum(A[i,i]).
    • transpose

      public DoubleMatrix2D transpose(DoubleMatrix2D A)
      Constructs and returns a new view which is the transposition of the given matrix A. Equivalent to A.viewDice(). This is a zero-copy transposition, taking O(1), i.e. constant time. The returned view is backed by this matrix, so changes in the returned view are reflected in this matrix, and vice-versa. Use idioms like result = transpose(A).copy() to generate an independent matrix.

      Example:

      2 x 3 matrix:
      1, 2, 3
      4, 5, 6
      transpose ==> 3 x 2 matrix:
      1, 4
      2, 5
      3, 6
      transpose ==> 2 x 3 matrix:
      1, 2, 3
      4, 5, 6
      Returns:
      a new transposed view.
    • trapezoidalLower

      protected DoubleMatrix2D trapezoidalLower(DoubleMatrix2D A)
      Modifies the matrix to be a lower trapezoidal matrix.
      Returns:
      A (for convenience only).
      See Also:
      • invalid reference
        #triangulateLower(DoubleMatrix2D)
    • xmultOuter

      private DoubleMatrix2D xmultOuter(DoubleMatrix1D x, DoubleMatrix1D y)
      Outer product of two vectors; Returns a matrix with A[i,j] = x[i] * y[j].
      Parameters:
      x - the first source vector.
      y - the second source vector.
      Returns:
      the outer product A.
    • xpowSlow

      private DoubleMatrix2D xpowSlow(DoubleMatrix2D A, int k)
      Linear algebraic matrix power; B = Ak invalid input: '<'==> B = A*A*...*A.
      Parameters:
      A - the source matrix; must be square.
      k - the exponent, can be any number.
      Returns:
      a new result matrix.
      Throws:
      IllegalArgumentException - if !Testing.isSquare(A).