Null.java

/*
 * $Id: Null.java,v 1.24 2008/03/15 00:23:40 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.matrix.IntMatrix;
import org.mklab.nfc.scalar.DoubleNumber;
import org.mklab.nfc.scalar.DoubleNumberUtil;


/**
 * 零空間の基底からなる行列を求めるクラスです。
 * 
 * <p>Null space basis
 * 
 * @author koga
 * @version $Revision: 1.24 $
 * @see org.mklab.tool.matrix.Orth
 */
public class Null {

  /**
   * <code>A</code> の零空間の直交基底を返します。 <code>Q</code> の列ベクトルが <code>A</code> の零空間を張り、 <code>Q</code> の列の数は、 <code>(Rows(A) - rank(A))</code> と等しい。
   * 
   * @param A 対象となる行列
   * @return 零空間の基底からなる行列 (null space basis)
   */
  public static DoubleMatrix nullMatrix(DoubleMatrix A) {
    double tolerance = DoubleNumberUtil.EPS * A.frobNorm().doubleValue();
    return nullMatrix(A, tolerance);
  }

  /**
   * <code>tolerance</code> は、零空間の次元を決めるために使われる。
   * 
   * @param A 対象となる行列
   * @param tolerance 許容誤差
   * @return 零空間の基底からなる行列 (null space basis)
   */
  public static DoubleMatrix nullMatrix(DoubleMatrix A, double tolerance) {
    QRDecomposition<DoubleNumber,DoubleMatrix> tmp = A.conjugateTranspose().qrDecomposeWithPermutation();
    DoubleMatrix Q = tmp.getQ();
    DoubleMatrix R = tmp.getR();
    // DoubleMatrix P = tmp[2];

    IntMatrix idx;
    if (Math.min(A.getColumnSize(), A.getRowSize()) > 1) {
      idx = R.absElementWise().sumRowWise().conjugateTranspose().compareElementWise(".<=", tolerance).find(); //$NON-NLS-1$
    } else {
      DoubleNumber unit = A.getElement(1, 1);
      idx = unit.createGrid(new DoubleNumber[] {R.getElement(1).abs()}).compareElementWise(".<=", tolerance).find(); //$NON-NLS-1$
    }

    if (idx.length() > 0) {
      Q = Q.getSubMatrix(1, Q.getRowSize(), idx);
    } else {
      Q = A.createZero(0, 0);
    }
    return Q;
  }

}