Orth.java

/*
 * $Id: Orth.java,v 1.20 2008/03/03 15:18:35 koga Exp $
 * 
 * Copyright (C) 2004 Koga Laboratory. All rights reserved.
 */
package org.mklab.tool.matrix;

import org.mklab.nfc.eig.QRDecomposition;
import org.mklab.nfc.matrix.DoubleMatrix;
import org.mklab.nfc.scalar.DoubleNumber;
import org.mklab.nfc.scalar.DoubleNumberUtil;


/**
 * 直交基底からなる行列を求めるクラスです。
 * 
 * @author koga
 * @version $Revision: 1.20 $
 * @see org.mklab.tool.matrix.Null
 */
public class Orth {

  /**
   * <code>A</code> の像の直交基底を返します。
   * 
   * <code>A</code> の像の直交基底の列の数は、 <code>rank(A)</code> と等しいです。
   * 
   * @param A 対象となる行列
   * @return 直交基底からなる行列 (orthogonal basis)
   */
  public static DoubleMatrix orth(DoubleMatrix A) {
    final double tolerance = DoubleNumberUtil.EPS * A.frobNorm().doubleValue();
    return orth(A, tolerance);
  }

  /**
   * <code>A</code> の像の直交基底を返します。
   * 
   * <code>A</code> の像の直交基底の列の数は、 <code>rank(A)</code> と等しいです。
   * 
   * <p><code>tolerance</code> は <code>A</code> のランクを決めるために使われます。
   * 
   * @param A 対象となる行列
   * @param tolerance 許容誤差
   * @return 直交基底からなる行列 (orthogonal basis)
   */
  public static DoubleMatrix orth(DoubleMatrix A, double tolerance) {
    final QRDecomposition<DoubleNumber,DoubleMatrix> qrp = A.qrDecomposeWithPermutation();
    DoubleMatrix Q = qrp.getQ();
    final DoubleMatrix R = qrp.getR();

    int rk = R.diagonalToVector().absElementWise().compareElementWise(".>", tolerance).find().length(); //$NON-NLS-1$

    if (rk > 0) {
      Q = Q.getColumnVectors(1, rk);
      Q = Q.unaryMinus();
      Q.setColumnVector(rk, Q.getColumnVector(rk).unaryMinus());
      return Q;
    }

    return A.createZero(0, 0);
  }
}