Mark Bishop

Analytical Chemistry and

Signal Acquisition/Analysis


/**
 * Copyright 2009 Mark Bishop
 *  This program is free software: you can redistribute it and/or modify
    it under the terms of the GNU General Public License as published by
    the Free Software Foundation, either version 3 of the License, or
    (at your option) any later version. 
 
    This program is distributed in the hope that it will be useful,
    but WITHOUT ANY WARRANTY; without even the implied warranty of
    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
    GNU General Public License for more details: http://www.gnu.org/licenses
 
    The author makes no warranty for the accuracy, completeness, safety, or usefulness
    of any information provided and does not represent that its use would not infringe
    privately owned right.
 */
 
import java.util.ArrayList;
import java.util.List;
 
public class CpuQR {
 
	/**
	 * Implements Golub/Van Loan Algorithm 5.1.1
	 * 
	 * @param vector
	 *            x[n]
	 * @return vBeta[n + 1] where: vBeta[0:n-1] contains the essential part of
	 *         the Householder vector (v[0] = 1) and vBeta[n] =
	 *         2/((vTranspose)v)
	 */
	public static float[] House(float[] x) {
		int m = x.length;
		float beta = 0;
		float[] vBeta = new float[m + 1];
 
		// scale to prevent overflow
		float N1 = CpuMatrixOps.VectorInfinityNorm(x);
		if (N1 != 0) {
			for (int i = 0; i < m; i++) {
				x[i] = x[i] / N1;
			}
		}
 
		float u = 0;
		float sigma = 0;
 
		for (int i = 1; i < m; i++) {
			sigma = (float) (sigma + Math.pow(x[i], 2));
			vBeta[i] = x[i];
		}
 
		if (sigma == 0) {
			beta = 0;
		} else {
			u = (float) Math.pow(((Math.pow(x[0], 2)) + sigma), 0.5);
			if (x[0] <= 0) {
				vBeta[0] = x[0] - u;
			} else {
				vBeta[0] = -sigma / (x[0] + u);
			}
			beta = (float) (2 * Math.pow(vBeta[0], 2) / (sigma + Math.pow(
					vBeta[0], 2)));
			for (int i = 1; i < m; i++) {
				vBeta[i] = vBeta[i] / vBeta[0];
			}
		}
		vBeta[m] = beta;
		vBeta[0] = 1;
 
		return vBeta;
	}
 
	/**
	 * Implements Golub/Van Loan 5.2.1
	 * 
	 * Overwrites source_matrix A with R.
	 * 
	 * This method is provided for illustrative purposes. It explicitly, and
	 * inefficiently, produces the matrices Q and R as complete and separate
	 * matrices.
	 * 
	 * @param A
	 *            A[m,n]: to be overwritten with R[m,n]
	 * @return Q[m,m] such that QR=A and Q is orthogonal
	 */
	public static float[][] HouseholderQR(float[][] A) {
		int m = A.length;
		int n = A[0].length;
		float[][] Q = CpuMatrixOps.Eye(m);
		float[] vector = new float[m];
		float[] vHouse = new float[m];
		float[] beta = new float[n];
 
		for (int j = 0; j < n; j++) {
			vHouse = new float[m];
			vector = new float[m - j];
			for (int i = 0; i < vector.length; i++) {
				vector[i] = A[j + i][j];
			}
			vector = House(vector);
			// copy vector into v, with leading zeros up to the jth row
			for (int i = 0; i < vector.length - 1; i++) {
				vHouse[j + i] = vector[i];
			}
			beta[j] = vector[vector.length - 1];
 
			float[] tmp1 = CpuMatrixOps.Multiply(CpuMatrixOps.Transpose(A),
					vHouse);
			float[][] omegaT = CpuMatrixOps.VectorOuterProduct(vHouse,
					CpuMatrixOps.Scale(beta[j], tmp1));
 
			for (int i = j; i < m; i++) {
				for (int k = j; k < n; k++) {
					A[i][k] = A[i][k] - omegaT[i][k];
				}
			}
			float[][] omega = CpuMatrixOps.Scale(
					beta[j],
					CpuMatrixOps.VectorOuterProduct(
							CpuMatrixOps.Multiply(Q, vHouse), vHouse));
 
			for (int i = 0; i < m; i++) {
				for (int k = 0; k < m; k++) {
					Q[i][k] = Q[i][k] - omega[i][k];
				}
			}
		}
		return Q;
	}
 
	// need Householder Compact QR = HouseholderCompactQR
 
	/**
	 * Implements Golub/ Van Loan 5.4.1
	 * 
	 * Overwrites source_matrix A with a compact form QR.
	 * 
	 * @param A
	 *            A[m,n] will be overwritten by QR
	 * @return an object list that contains: The matrix rank (int), the
	 *         permutation vector(int[]), and beta(float[]). Beta is a vector of
	 *         factors such that Beta[j] = 2/(v[j]Transpose *v[j]) and
	 *         Beta[length-a] = 0. Beta is not a required element of a QR Pivot
	 *         factorization because it can be re-calculated as needed from the
	 *         compact QR, but consuming this vector may eliminate the overhead
	 *         of re-computing the values in subsequent functions.
	 */
	public static List<Object> HouseholderCompactQRPivot(float[][] A) {
		int m = A.length;
		int n = A[0].length;
		int tIndex = 0;
		int rankIndex;
		int rank = 0;
		int[] piv = new int[n];
		float tau = (float) 0.000000000000; // Minimum norm threshold
		float maxMag = tau;
		float[] beta = new float[n];
		float[] mag = new float[n];
 
		for (int i = 0; i < n; i++) {
			piv[i] = i;
		}
 
		// A vector of column magnitudes = mag[j] = A(1:m,j)Transpose * A(1:m,j)
		// =
		// 2-norm**2
		for (int j = 0; j < n; j++) {
			for (int i = 0; i < m; i++) {
				mag[j] = (float) (mag[j] + (Math.pow(A[i][j], 2)));
			}
		}
 
		// find column with highest magnitude and store that 2-norm**2 in maxMag
		for (int i = 0; i < n; i++) {
			if ((mag[i]) > maxMag) {
				maxMag = mag[i];
				tIndex = i;
			}
		}
 
		while (maxMag > tau) {
			rank = rank + 1;
			rankIndex = rank - 1;
			// column swap
			if (tIndex + 1 != rank) {
				// for the vector of magnitudes
				float[] ColTemp = new float[m];
				float cTemp = 0;
				cTemp = mag[tIndex];
				mag[tIndex] = mag[rankIndex];
				mag[rankIndex] = cTemp;
				// for the A matrix
				for (int i = 0; i < m; i++) {
					ColTemp[i] = A[i][rankIndex];
					A[i][rankIndex] = A[i][tIndex];
					A[i][tIndex] = ColTemp[i];
				}
				// for the permutation vector
				int temp = 0;
				temp = piv[tIndex];
				piv[tIndex] = piv[rankIndex];
				piv[rankIndex] = temp;
			}
 
			// Extract the column with index = rankIndex to a float[] for
			// submission to House.
			float[] vector = new float[m - rankIndex];
			for (int i = 0; i < vector.length; i++) {
				vector[i] = A[rankIndex + i][rankIndex];
			}
 
			vector = House(vector);
 
			// Convert the result to a column vector with leading zeros up to
			// the row at rankIndex.
			float[] vHouse = new float[m];
			for (int i = 0; i < vector.length - 1; i++) {
				vHouse[rankIndex + i] = vector[i];
			}
 
			// Accumulate R
			beta[rankIndex] = vector[vector.length - 1];
			float[] tmp1 = CpuMatrixOps.Multiply(CpuMatrixOps.Transpose(A),
					vHouse);
			float[][] omegaT = CpuMatrixOps.VectorOuterProduct(vHouse,
					CpuMatrixOps.Scale(beta[rankIndex], tmp1));
 
			for (int i = rankIndex; i < m; i++) {
				for (int j = rankIndex; j < n; j++) {
					A[i][j] = A[i][j] - omegaT[i][j];
				}
			}
 
			// Place the essential part of the Householder vectors below the
			// principal diagonal
			for (int i = rankIndex + 1; i < m; i++) {
				A[i][rankIndex] = vector[i - rankIndex];
			}
 
			// Re-evaluate mag[]
			for (int i = rankIndex + 1; i < n; i++) {
				mag[i] = mag[i] - (A[rankIndex][i] * A[rankIndex][i]);
			}
			if (rank < n) {
				maxMag = tau;
				for (int i = rankIndex + 1; i < n; i++) {
					if (mag[i] > maxMag) {
						maxMag = mag[i];
						tIndex = i;
					}
				}
			} else {
				maxMag = 0;
			}
		}
 
		List<Object> oPBeta = new ArrayList<Object>();
		oPBeta.add(rank);
		oPBeta.add(piv);
		oPBeta.add(beta);
 
		return oPBeta;
	}
 
	/**
	 * Parses compact QR into Q and R and may take Beta, a vector where Beta[j]
	 * =v[j]Transpose*v[j] and Beta[length-1] = 0
	 * 
	 * @param Encoded_matrix
	 * @param beta
	 *            optional pre-calculated Beta vector
	 * @return Q and overwrites QR with R
	 */
	public static float[][] ParseHouseholderQR(float[][] QR, float[] beta) {
		int m = QR.length;
		int n = QR[0].length;
		float[][] Q = CpuMatrixOps.Eye(m);
 
		for (int j = n - 1; j >= 0; j += -1) {
			float[] v = new float[m];
			v[j] = 1;
			for (int i = j + 1; i < m; i++) {
				v[i] = QR[i][j];
			}
 
			float b = 0;
			if (beta != null) {
				b = beta[j];
			} else {
				if (j < n - 1) {
					b = 2 / CpuMatrixOps.VectorInnerProduct(v, v);
				}
			}
 
			float[] omega = CpuMatrixOps.Multiply(CpuMatrixOps.Transpose(Q), v);
			for (int i = 0; i < m; i++) {
				for (int k = 0; k < m; k++) {
					Q[i][k] = Q[i][k] - b * omega[k] * v[i];
				}
			}
		}
		for (int j = 0; j < n; j++) {
			for (int i = j + 1; i < m; i++) {
				QR[i][j] = 0;
			}
		}
		return Q;
	}
 
	/**
	 * Generates a permutation matrix from a permutation vector. This method is
	 * provided for illustrative purposes.Efficient methods do not store/use
	 * permutation matrices. Instead, the permutation matrix is stored in
	 * compact vector form. It's elements are column indices and provide a swap
	 * order. For example, the permutation vector {1, 0} represents a
	 * permutation matrix that consists of I[2] with the columns swapped, i.e.:
	 * {{0, 1}, {1, 0}}.
	 * 
	 * @param P
	 *            Permutation vector
	 * @return Permutation matrix
	 */
	public static float[][] ParseColumnPermutation(int[] P) {
		int n = P.length;
		float[][] Q = new float[n][n];
		float[][] identity = CpuMatrixOps.Eye(n);
		for (int i = 0; i < n; i++) {
			for (int j = 0; j < n; j++) {
				Q[i][j] = identity[i][P[j]];
			}
		}
		return Q;
	}
 
	/**
	 * Golub/Van Loan pp 259, 239. Given Q, R, P, rank, b; solve for x: Ax = b
	 * where AP = QR
	 * 
	 * @param Q
	 * @param R
	 * @param P
	 * @param rank
	 * @param b
	 * @return x: Ax = b
	 */
	public static float[] QRPivotSolveForX(float[][] Q, float[][] R, int[] P,
			int rank, float[] b) {
		int m = Q.length;
		int n = R[0].length;
 
		// make QTranspose*b = QTb. A thin Q is used here.
		float[] QTb = new float[rank]; // QTb=c
		for (int i = 0; i < rank; i++) {
			for (int j = 0; j < m; j++) {
				QTb[i] = Q[j][i] * b[j] + QTb[i]; // Q(j,i)=QT(i,j)
			}
		}
		// By back calculation QTAPz=QTb and QTAP is upper triangular. The first
		// rank X rank rows and cols of R = QTAP=r11
		float[] z = new float[n];
		z[rank - 1] = QTb[rank - 1] / R[rank - 1][rank - 1];
		for (int i = rank - 2; i >= 0; i += -1) {
			for (int j = rank - 1; j >= i + 1; j += -1) {
				z[i] = z[i] - R[i][j] * z[j];
			}
			z[i] = z[i] + QTb[i];
			z[i] = z[i] / R[i][i];
		}
 
		// Apply permutation vector
		float[] tempb = new float[n];
		for (int i = 0; i < n; i++) {
			tempb[P[i]] = z[i];
		}
		return tempb;
	}
 
	/**
	 * Golub/Van Loan pp 259, 239. Given: a compact QR, Beta, P, rank, b; solve
	 * for x: Ax = b where AP = QR. If Beta = null, then Beta will be
	 * calculated. The permutation vector P may be (QR without pivoting).
	 * 
	 * @param QR
	 * @param beta
	 * @param P
	 * @param rank
	 * @param b
	 * @return x: Ax = b
	 */
	public static float[] HouseholderCompactQRSolveForX(float[][] QR,
			float[] beta, int[] P, int rank, float[] b) {
		int m = QR.length;
		int n = QR[0].length;
		float[] QTb = Utilities.CloneVector(b);
 
		// Hb=(I-BetavvT)b=A-v(BetabTv)T Golub/Van Loan p 211
		for (int j = 0; j < rank; j++) {
			float[] v = new float[m];
			v[j] = 1;
			for (int i = j + 1; i < m; i++) {
				v[i] = QR[i][j];
			}
 
			float omega = CpuMatrixOps.VectorInnerProduct(QTb, v);
 
			float B = 0;
			if (beta != null) {
				B = beta[j];
			} else {
				if (j < n - 1) {
					B = 2 / CpuMatrixOps.VectorInnerProduct(v, v);
				}
			}
			for (int i = 0; i < m; i++) {
				QTb[i] = QTb[i] - B * omega * v[i];
			}
		}
 
		// QTAPx=QTb and QTAP is the upper triangular m=n=rank = rank super
		// matrix of R
		// (r11).
		float[] z = new float[n];
		z[rank - 1] = QTb[rank - 1] / QR[rank - 1][rank - 1];
		for (int i = rank - 2; i >= 0; i += -1) {
			for (int j = rank - 1; j >= i + 1; j += -1) {
				z[i] = z[i] - QR[i][j] * z[j];
			}
			z[i] = z[i] + QTb[i];
			z[i] = z[i] / QR[i][i];
		}
 
		// Apply permutation vector
		if (P != null) {
			float[] tempb = new float[n];
			for (int i = 0; i < n; i++) {
				tempb[P[i]] = z[i];
 
			}
			return tempb;
		}
 
		return z;
	}
 
	/**
	 * QR Golub/Van Loan 5.2.1 Form encoded QR without explicit formation of
	 * Householder matrices
	 * 
	 * @param A
	 * @return Beta and overwrites A with compact QR
	 */
	public static float[] HouseholderQREncoded(float[][] A) {
		int m = A.length;
		int n = A[0].length;
		float[] beta = new float[n];
		for (int j = 0; j < n; j++) {
			float[] vector = new float[m - j];
			for (int i = 0; i < vector.length; i++) {
				vector[i] = A[j + i][j];
			}
			vector = House(vector);
 
			// Convert the result to a column vector with leading zeros up to
			// the row at rankIndex.
			float[] vHouse = new float[m];
			for (int i = 0; i < vector.length - 1; i++) {
				vHouse[j + i] = vector[i];
			}
 
			// Accumulate R
			beta[j] = vector[vector.length - 1];
			float[] tmp1 = CpuMatrixOps.Multiply(CpuMatrixOps.Transpose(A),
					vHouse);
			float[][] omegaT = CpuMatrixOps.VectorOuterProduct(vHouse,
					CpuMatrixOps.Scale(beta[j], tmp1));
 
			for (int i = j; i < m; i++) {
				for (int k = j; k < n; k++) {
					A[i][k] = A[i][k] - omegaT[i][k];
				}
			}
			// populate the lower triangular matrix holding the Householder
			// vectors
			for (int i = j + 1; i < m; i++) {
				A[i][j] = vector[i - j];
			}
		}
 
		return beta;
	}
}
 
 
ID may be exposed