matrix.java

一个数据挖掘软件ALPHAMINERR的整个过程的JAVA版源代码

  
  /**
   * Returns true if the matrix is symmetric.
   *
   * @return boolean true if matrix is symmetric.
   */
  public boolean isSymmetric() {

    int nr = m_Elements.length, nc = m_Elements[0].length;
    if (nr != nc)
      return false;

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

  /**
   * Returns the multiplication of two matrices
   *
   * @param b the multiplication matrix
   * @return the product matrix
   */
  public final Matrix multiply(Matrix b) {
   
    int nr = m_Elements.length, nc = m_Elements[0].length;
    int bnr = b.m_Elements.length, bnc = b.m_Elements[0].length;
    Matrix c = new Matrix(nr, bnc);

    for(int i = 0; i < nr; i++) {
      for(int j = 0; j< bnc; j++) {
	for(int k = 0; k < nc; k++) {
	  c.m_Elements[i][j] += m_Elements[i][k] * b.m_Elements[k][j];
	}
      }
    }

    return c;
  }

  /**
   * Performs a (ridged) linear regression.
   *
   * @param y the dependent variable vector
   * @param ridge the ridge parameter
   * @return the coefficients 
   * @exception IllegalArgumentException if not successful
   */
  public final double[] regression(Matrix y, double ridge) {

    if (y.numColumns() > 1) {
      throw new IllegalArgumentException("Only one dependent variable allowed");
    }
    int nc = m_Elements[0].length;
    double[] b = new double[nc];
    Matrix xt = this.transpose();

    boolean success = true;

    do {
      Matrix ss = xt.multiply(this);
      
      // Set ridge regression adjustment
      for (int i = 0; i < nc; i++) {
	ss.setElement(i, i, ss.getElement(i, i) + ridge);
      }

      // Carry out the regression
      Matrix bb = xt.multiply(y);
      for(int i = 0; i < nc; i++) {
	b[i] = bb.m_Elements[i][0];
      }
      try {
	ss.solve(b);
	success = true;
      } catch (Exception ex) {
	ridge *= 10;
	success = false;
      }
    } while (!success);

    return b;
  }

  /**
   * Performs a weighted (ridged) linear regression. 
   *
   * @param y the dependent variable vector
   * @param w the array of data point weights
   * @param ridge the ridge parameter
   * @return the coefficients 
   * @exception IllegalArgumentException if the wrong number of weights were
   * provided.
   */
  public final double[] regression(Matrix y, double [] w, double ridge) {

    if (w.length != numRows()) {
      throw new IllegalArgumentException("Incorrect number of weights provided");
    }
    Matrix weightedThis = new Matrix(numRows(), numColumns());
    Matrix weightedDep = new Matrix(numRows(), 1);
    for (int i = 0; i < w.length; i++) {
      double sqrt_weight = Math.sqrt(w[i]);
      for (int j = 0; j < numColumns(); j++) {
	weightedThis.setElement(i, j, getElement(i, j) * sqrt_weight);
      }
      weightedDep.setElement(i, 0, y.getElement(i, 0) * sqrt_weight);
    }
    return weightedThis.regression(weightedDep, ridge);
  }

  /**
   * Returns the L part of the matrix.
   * This does only make sense after LU decomposition.
   *
   * @return matrix with the L part of the matrix; 
   */
  public Matrix getL() throws Exception {

    int nr = m_Elements.length;    // num of rows
    int nc = m_Elements[0].length; // num of columns
    double[][] ld = new double[nr][nc];
    
    for (int i = 0; i < nr; i++) {
      for (int j = 0; (j < i) && (j < nc); j++) {
	ld[i][j] = m_Elements[i][j];
      }
      if (i < nc) ld[i][i] = 1;
    }
    Matrix l = new Matrix(ld);
    return l;
  }

  /**
   * Returns the U part of the matrix.
   * This does only make sense after LU decomposition.
   *
   * @return matrix with the U part of a matrix; 
   */
  public Matrix getU() throws Exception {

    int nr = m_Elements.length;    // num of rows
    int nc = m_Elements[0].length; // num of columns
    double[][] ud = new double[nr][nc];
    
    for (int i = 0; i < nr; i++) {
      for (int j = i; j < nc ; j++) {
	ud[i][j] = m_Elements[i][j];
      }
    }
    Matrix u = new Matrix(ud);
    return u;
  }
  
  /**
   * Performs a LUDecomposition on the matrix.
   * It changes the matrix into its LU decomposition using the
   * Crout algorithm
   *
   * @return the indices of the row permutation
   * @exception Exception if the matrix is singular 
   */
  public int [] LUDecomposition() throws Exception {
    
    int nr = m_Elements.length;    // num of rows
    int nc = m_Elements[0].length; // num of columns

    int [] piv = new int[nr];

    double [] factor = new double[nr];
    double dd;
    for (int row = 0; row < nr; row++) {
      double max = Math.abs(m_Elements[row][0]);
      for (int col = 1; col < nc; col++) {
	if ((dd = Math.abs(m_Elements[row][col])) > max)
	  max = dd;
      }
      if (max < 0.000000001) {
	throw new Exception("Matrix is singular!");
      }
      factor[row] = 1 / max;
    }

    for (int i = 1; i < nr; i++) piv[i] = i;

    for (int col = 0; col < nc; col++) {

      // compute beta i,j
      for (int row = 0; (row <= col) && (row < nr); row ++) {
	double sum = 0.0;

	for (int k = 0; k < row; k++) {
	  sum += m_Elements[row][k] * m_Elements[k][col];
	}
	m_Elements[row][col] = m_Elements[row][col] - sum; // beta i,j
;
      }

      // compute alpha i,j
      for (int row = col + 1; row < nr; row ++) {
	double sum = 0.0;
	for (int k = 0; k < col; k++) {

	  sum += m_Elements[row][k] * m_Elements[k][col];
	}
        // alpha i,j before division
	m_Elements[row][col] = (m_Elements[row][col] - sum);

      }
      
      // find row for division:see if any of the alphas are larger then b j,j
      double maxFactor = Math.abs(m_Elements[col][col]) * factor[col];
      int newrow = col;

      for (int ii = col + 1; ii < nr; ii++) {
	if ((Math.abs(m_Elements[ii][col]) * factor[ii]) > maxFactor) {
	  newrow = ii; // new row
	  maxFactor = Math.abs(m_Elements[ii][col]) * factor[ii];
	}  
      }
      if (newrow != col) {
        // swap the rows
	for (int kk = 0; kk < nc; kk++) {
	  double hh = m_Elements[col][kk];
	  m_Elements[col][kk] = m_Elements[newrow][kk];
	  m_Elements[newrow][kk] = hh;
	}
	// remember pivoting
	int help = piv[col]; piv[col] = piv[newrow]; piv[newrow] = help;
	double hh  = factor[col]; factor[col] = factor[newrow]; factor[newrow] = help;
      }

      if (m_Elements[col][col] == 0.0) {
	throw new Exception("Matrix is singular");
      }

      for (int row = col + 1; row < nr; row ++) {
	// divide all 
	m_Elements[row][col] = m_Elements[row][col] / 
	  m_Elements[col][col];
      }
   }

    return piv;
  }
  
  /**
   * Solve A*X = B using backward substitution.
   * A is current object (this). Note that this matrix will be changed! 
   * B parameter bb
   * X returned in parameter bb
   *
   * @param bb first vektor B in above equation then X in same equation.
   * @exception matrix is singulaer
   */
  public void solve(double [] bb) throws Exception {

    int nr = m_Elements.length;
    int nc = m_Elements[0].length;

    double [] b = new double[bb.length];
    double [] x = new double[bb.length];
    double [] y = new double[bb.length];

    int [] piv = this.LUDecomposition();

    // change b according to pivot vector
    for (int i = 0; i < piv.length; i++) {
      b[i] = bb[piv[i]];
    }
    
    y[0] = b[0];

    for (int row = 1; row < nr; row++) {
      double sum = 0.0;
      for (int col = 0; col < row; col++) {

	sum += m_Elements[row][col] * y[col];
      }
      y[row] = b[row] - sum;
    }

    x[nc - 1] = y[nc - 1] / m_Elements[nc - 1][nc - 1];
    for (int row = nc - 2; row >= 0; row--) {
      double sum = 0.0;
      for (int col = row + 1; col < nc; col++) {

	sum += m_Elements[row][col] * x[col];
      }
      x[row] = (y[row] - sum) / m_Elements[row][row];
    }

    // change x according to pivot vector
    for (int i = 0; i < piv.length; i++) {

      //bb[piv[i]] = x[i];
      bb[i] = x[i];
    }
  }

  /**
   * Performs Eigenvalue Decomposition using Householder QR Factorization
   *
   * This function is adapted from the CERN Jet Java libraries, for it 
   * the following copyright applies (see also, text on top of file)
   * Copyright (C) 1999 CERN - European Organization for Nuclear Research.
   *
   * Matrix must be symmetrical.
   * Eigenvectors are return in parameter V, as columns of the 2D array.
   * Eigenvalues are returned in parameter d.
   *
   *
   * @param V double array in which the eigenvectors are returned 
   * @param d array in which the eigenvalues are returned
   * @exception if matrix is not symmetric
   */

  
  public void eigenvalueDecomposition(double [][] V, double [] d) throws Exception {

    if (!this.isSymmetric())
      throw new Exception("EigenvalueDecomposition: Matrix must be symmetric.");
    
    int n = this.numRows();
    double[]  e = new double [n]; //todo, don't know what this e is really for!
    
    for (int i = 0; i < n; i++) {
      for (int j = 0; j < n; j++) {
	V[i][j] = m_Elements[i][j];
      }
    }
      
    // Tridiagonalize using householder reduction
    tred2(V, d, e, n);
      
    // Diagonalize.
    tql2(V, d, e, n);

    // Matrix v = new Matrix (V);
    // testEigen(v, d);
  } 


  /**
   * Symmetric Householder reduction to tridiagonal form.
   *
   * This function is adapted from the CERN Jet Java libraries, for it 
   * the following copyright applies (see also, text on top of file)
   * Copyright (C) 1999 CERN - European Organization for Nuclear Research.
   *
   * @param V containing  a copy of the values of the matrix
   * @param d 
   * @param e
   * @param n size of matrix
   * @exception if reduction doesn't work
   */
  private void tred2 (double [][] V, double [] d, double [] e, int n) throws Exception {
    //  This is derived from the Algol procedures tred2 by
    //  Bowdler, Martin, Reinsch, and Wilkinson, Handbook for
    //  Auto. Comp., Vol.ii-Linear Algebra, and the corresponding
    //  Fortran subroutine in EISPACK.
    
    for (int j = 0; j < n; j++) {
      d[j] = V[n-1][j];
    }
    
    // Householder reduction to tridiagonal form.

    // Matrix v = new Matrix(V);
    // System.out.println("before household - Matrix V \n" + v);
    for (int i = n-1; i > 0; i--) {
   
      // Matrix v = new Matrix(V);
      // System.out.println("Matrix V \n" + v);

      // Scale to avoid under/overflow.
      double scale = 0.0;
      double h = 0.0;
      for (int k = 0; k < i; k++) {
	scale = scale + Math.abs(d[k]);
      }
      if (scale == 0.0) {
	e[i] = d[i-1];
	for (int j = 0; j < i; j++) {
	  d[j] = V[i-1][j];
	  V[i][j] = 0.0;
	  V[j][i] = 0.0;
	}