matrix.java

MacroWeka扩展了著名数据挖掘工具weka

    try {
      for (int i = 0; i < r.length; i++) {
        for (int j = j0; j <= j1; j++) {
          A[r[i]][j] = X.get(i,j-j0);
        }
      }
    } catch(ArrayIndexOutOfBoundsException e) {
      throw new ArrayIndexOutOfBoundsException("Submatrix indices");
    }
  }

  /** 
   * Set a submatrix.
   * @param i0   Initial row index
   * @param i1   Final row index
   * @param c    Array of column indices.
   * @param X    A(i0:i1,c(:))
   * @exception  ArrayIndexOutOfBoundsException Submatrix indices
   */
  public void setMatrix(int i0, int i1, int[] c, Matrix X) {
    try {
      for (int i = i0; i <= i1; i++) {
        for (int j = 0; j < c.length; j++) {
          A[i][c[j]] = X.get(i-i0,j);
        }
      }
    } catch(ArrayIndexOutOfBoundsException e) {
      throw new ArrayIndexOutOfBoundsException("Submatrix indices");
    }
  }
  
  /**
   * Returns true if the matrix is symmetric.
   *
   * @return boolean true if matrix is symmetric.
   * @author FracPete, taken from old weka.core.Matrix class
   */
  public boolean isSymmetric() {
    int nr = A.length, nc = A[0].length;
    if (nr != nc)
      return false;

    for (int i = 0; i < nc; i++) {
      for (int j = 0; j < i; j++) {
        if (A[i][j] != A[j][i])
          return false;
      }
    }
    return true;
  }

  /**
   * returns whether the matrix is a square matrix or not.
   *
   * @return true if the matrix is a square matrix
   * @author FracPete
   */
  public boolean isSquare() {
    return (getRowDimension() == getColumnDimension());
  }

  /** 
   * Matrix transpose.
   * @return    A'
   */
  public Matrix transpose() {
    Matrix X = new Matrix(n,m);
    double[][] C = X.getArray();
    for (int i = 0; i < m; i++) {
      for (int j = 0; j < n; j++) {
        C[j][i] = A[i][j];
      }
    }
    return X;
  }

  /** 
   * One norm
   * @return    maximum column sum.
   */
  public double norm1() {
    double f = 0;
    for (int j = 0; j < n; j++) {
      double s = 0;
      for (int i = 0; i < m; i++) {
        s += Math.abs(A[i][j]);
      }
      f = Math.max(f,s);
    }
    return f;
  }

  /** 
   * Two norm
   * @return    maximum singular value.
   */
  public double norm2() {
    return (new SingularValueDecomposition(this).norm2());
  }

  /** 
   * Infinity norm
   * @return    maximum row sum.
   */
  public double normInf() {
    double f = 0;
    for (int i = 0; i < m; i++) {
      double s = 0;
      for (int j = 0; j < n; j++) {
        s += Math.abs(A[i][j]);
      }
      f = Math.max(f,s);
    }
    return f;
  }

  /** 
   * Frobenius norm
   * @return    sqrt of sum of squares of all elements.
   */
  public double normF() {
    double f = 0;
    for (int i = 0; i < m; i++) {
      for (int j = 0; j < n; j++) {
        f = Maths.hypot(f,A[i][j]);
      }
    }
    return f;
  }

  /**  
   * Unary minus
   * @return    -A
   */
  public Matrix uminus() {
    Matrix X = new Matrix(m,n);
    double[][] C = X.getArray();
    for (int i = 0; i < m; i++) {
      for (int j = 0; j < n; j++) {
        C[i][j] = -A[i][j];
      }
    }
    return X;
  }

  /** 
   * C = A + B
   * @param B    another matrix
   * @return     A + B
   */
  public Matrix plus(Matrix B) {
    checkMatrixDimensions(B);
    Matrix X = new Matrix(m,n);
    double[][] C = X.getArray();
    for (int i = 0; i < m; i++) {
      for (int j = 0; j < n; j++) {
        C[i][j] = A[i][j] + B.A[i][j];
      }
    }
    return X;
  }

  /** 
   * A = A + B
   * @param B    another matrix
   * @return     A + B
   */
  public Matrix plusEquals(Matrix B) {
    checkMatrixDimensions(B);
    for (int i = 0; i < m; i++) {
      for (int j = 0; j < n; j++) {
        A[i][j] = A[i][j] + B.A[i][j];
      }
    }
    return this;
  }

  /** 
   * C = A - B
   * @param B    another matrix
   * @return     A - B
   */
  public Matrix minus(Matrix B) {
    checkMatrixDimensions(B);
    Matrix X = new Matrix(m,n);
    double[][] C = X.getArray();
    for (int i = 0; i < m; i++) {
      for (int j = 0; j < n; j++) {
        C[i][j] = A[i][j] - B.A[i][j];
      }
    }
    return X;
  }

  /** 
   * A = A - B
   * @param B    another matrix
   * @return     A - B
   */
  public Matrix minusEquals(Matrix B) {
    checkMatrixDimensions(B);
    for (int i = 0; i < m; i++) {
      for (int j = 0; j < n; j++) {
        A[i][j] = A[i][j] - B.A[i][j];
      }
    }
    return this;
  }

  /** 
   * Element-by-element multiplication, C = A.*B
   * @param B    another matrix
   * @return     A.*B
   */
  public Matrix arrayTimes(Matrix B) {
    checkMatrixDimensions(B);
    Matrix X = new Matrix(m,n);
    double[][] C = X.getArray();
    for (int i = 0; i < m; i++) {
      for (int j = 0; j < n; j++) {
        C[i][j] = A[i][j] * B.A[i][j];
      }
    }
    return X;
  }

  /** 
   * Element-by-element multiplication in place, A = A.*B
   * @param B    another matrix
   * @return     A.*B
   */
  public Matrix arrayTimesEquals(Matrix B) {
    checkMatrixDimensions(B);
    for (int i = 0; i < m; i++) {
      for (int j = 0; j < n; j++) {
        A[i][j] = A[i][j] * B.A[i][j];
      }
    }
    return this;
  }

  /** 
   * Element-by-element right division, C = A./B
   * @param B    another matrix
   * @return     A./B
   */
  public Matrix arrayRightDivide(Matrix B) {
    checkMatrixDimensions(B);
    Matrix X = new Matrix(m,n);
    double[][] C = X.getArray();
    for (int i = 0; i < m; i++) {
      for (int j = 0; j < n; j++) {
        C[i][j] = A[i][j] / B.A[i][j];
      }
    }
    return X;
  }

  /** 
   * Element-by-element right division in place, A = A./B
   * @param B    another matrix
   * @return     A./B
   */
  public Matrix arrayRightDivideEquals(Matrix B) {
    checkMatrixDimensions(B);
    for (int i = 0; i < m; i++) {
      for (int j = 0; j < n; j++) {
        A[i][j] = A[i][j] / B.A[i][j];
      }
    }
    return this;
  }

  /** 
   * Element-by-element left division, C = A.\B
   * @param B    another matrix
   * @return     A.\B
   */
  public Matrix arrayLeftDivide(Matrix B) {
    checkMatrixDimensions(B);
    Matrix X = new Matrix(m,n);
    double[][] C = X.getArray();
    for (int i = 0; i < m; i++) {
      for (int j = 0; j < n; j++) {
        C[i][j] = B.A[i][j] / A[i][j];
      }
    }
    return X;
  }

  /** 
   * Element-by-element left division in place, A = A.\B
   * @param B    another matrix
   * @return     A.\B
   */
  public Matrix arrayLeftDivideEquals(Matrix B) {
    checkMatrixDimensions(B);
    for (int i = 0; i < m; i++) {
      for (int j = 0; j < n; j++) {
        A[i][j] = B.A[i][j] / A[i][j];
      }
    }
    return this;
  }

  /** 
   * Multiply a matrix by a scalar, C = s*A
   * @param s    scalar
   * @return     s*A
   */
  public Matrix times(double s) {
    Matrix X = new Matrix(m,n);
    double[][] C = X.getArray();
    for (int i = 0; i < m; i++) {
      for (int j = 0; j < n; j++) {
        C[i][j] = s*A[i][j];
      }
    }
    return X;
  }

  /** 
   * Multiply a matrix by a scalar in place, A = s*A
   * @param s    scalar
   * @return     replace A by s*A
   */
  public Matrix timesEquals(double s) {
    for (int i = 0; i < m; i++) {
      for (int j = 0; j < n; j++) {
        A[i][j] = s*A[i][j];
      }
    }
    return this;
  }

  /** 
   * Linear algebraic matrix multiplication, A * B
   * @param B    another matrix
   * @return     Matrix product, A * B
   * @exception  IllegalArgumentException Matrix inner dimensions must agree.
   */
  public Matrix times(Matrix B) {
    if (B.m != n) {
      throw new IllegalArgumentException("Matrix inner dimensions must agree.");
    }
    Matrix X = new Matrix(m,B.n);
    double[][] C = X.getArray();
    double[] Bcolj = new double[n];
    for (int j = 0; j < B.n; j++) {
      for (int k = 0; k < n; k++) {
        Bcolj[k] = B.A[k][j];
      }
      for (int i = 0; i < m; i++) {
        double[] Arowi = A[i];
        double s = 0;
        for (int k = 0; k < n; k++) {
          s += Arowi[k]*Bcolj[k];
        }
        C[i][j] = s;
      }
    }
    return X;
  }

  /** 
   * LU Decomposition
   * @return     LUDecomposition
   * @see LUDecomposition
   */
  public LUDecomposition lu() {
    return new LUDecomposition(this);
  }

  /** 
   * QR Decomposition
   * @return     QRDecomposition
   * @see QRDecomposition
   */
  public QRDecomposition qr() {
    return new QRDecomposition(this);
  }

  /** 
   * Cholesky Decomposition
   * @return     CholeskyDecomposition
   * @see CholeskyDecomposition
   */
  public CholeskyDecomposition chol() {
    return new CholeskyDecomposition(this);
  }

  /** 
   * Singular Value Decomposition
   * @return     SingularValueDecomposition
   * @see SingularValueDecomposition
   */
  public SingularValueDecomposition svd() {
    return new SingularValueDecomposition(this);
  }

  /** 
   * Eigenvalue Decomposition
   * @return     EigenvalueDecomposition
   * @see EigenvalueDecomposition
   */
  public EigenvalueDecomposition eig() {
    return new EigenvalueDecomposition(this);
  }

  /** 
   * Solve A*X = B
   * @param B    right hand side
   * @return     solution if A is square, least squares solution otherwise
   */
  public Matrix solve(Matrix B) {
    return (m == n ? (new LUDecomposition(this)).solve(B) :
        (new QRDecomposition(this)).solve(B));
  }

  /** 
   * Solve X*A = B, which is also A'*X' = B'
   * @param B    right hand side
   * @return     solution if A is square, least squares solution otherwise.
   */
  public Matrix solveTranspose(Matrix B) {
    return transpose().solve(B.transpose());
  }

  /** 
   * Matrix inverse or pseudoinverse
   * @return     inverse(A) if A is square, pseudoinverse otherwise.
   */
  public Matrix inverse() {
    return solve(identity(m,m));
  }

  /**
   * returns the square root of the matrix, i.e., X from the equation
   * X*X = A.<br/>
   * Steps in the Calculation (see <a href="http://www.mathworks.com/access/helpdesk/help/techdoc/ref/sqrtm.html" target="blank"><code>sqrtm</code></a> in Matlab):<br/>
   * <ol>
   *   <li>perform eigenvalue decomposition<br/>[V,D]=eig(A)</li>
   *   <li>take the square root of all elements in D (only the ones with 
   *       positive sign are considered for further computation)<br/>
   *       S=sqrt(D)</li>
   *   <li>calculate the root<br/>
   *       X=V*S/V, which can be also written as X=(V'\(V*S)')'</li>
   * </ol>
   * <p/>
   * <b>Note:</b> since this method uses other high-level methods, it generates
   * several instances of matrices. This can be problematic with large
   * matrices.
   * <p/>
   * Examples:
   * <ol>
   *   <li>
   *   <pre>
   *  X =
   *   5   -4    1    0    0
   *  -4    6   -4    1    0
   *   1   -4    6   -4    1
   *   0    1   -4    6   -4
   *   0    0    1   -4    5
   * 
   *  sqrt(X) =
   *   2   -1   -0   -0   -0 
   *  -1    2   -1    0   -0 
   *   0   -1    2   -1    0 
   *  -0    0   -1    2   -1 
   *  -0   -0   -0   -1    2 
   *  
   *  Matrix m = new Matrix(new double[][]{{5,-4,1,0,0},{-4,6,-4,1,0},{1,-4,6,-4,1},{0,1,-4,6,-4},{0,0,1,-4,5}});
   *   </pre>
   *   </li>
   *   <li>
   *   <pre>
   *  X =
   *   7   10
   *  15   22
   *  
   *  sqrt(X) =
   *  1.5667    1.7408
   *  2.6112    4.1779
   * 
   *  Matrix m = new Matrix(new double[][]{{7, 10},{15, 22}});
   *   </pre>
   *   </li>
   * </ol>
   *
   * @return    sqrt(A)
   * @author    FracPete
   */
  public Matrix sqrt() {
    EigenvalueDecomposition   evd;
    Matrix                    s;
    Matrix                    v;
    Matrix                    d;
    Matrix                    result;
    Matrix                    a;
    Matrix                    b;
    int                       i;
    int                       n;

    result = null;
    
    // eigenvalue decomp.
    // [V, D] = eig(A) with A = this
    evd = this.eig();
    v   = evd.getV();
    d   = evd.getD();

    // S = sqrt of cells of D
    s = new Matrix(d.getRowDimension(), d.getColumnDimension());