ptolemaeus.math.linear.householderreductions.ComplexHessenbergReduction Maven / Gradle / Ivy
/**
* MIT License
*
* Copyright (c) 2023 Lockheed Martin Corporation
*
* Permission is hereby granted, free of charge, to any person obtaining a copy of this software and associated
* documentation files (the "Software"), to deal in the Software without restriction, including without limitation the
* rights to use, copy, modify, merge, publish, distribute, sublicense, and/or sell copies of the Software, and to
* permit persons to whom the Software is furnished to do so, subject to the following conditions:
*
* The above copyright notice and this permission notice shall be included in all copies or substantial portions of the
* Software.
*
* THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE
* WARRANTIES OF MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE AUTHORS OR
* COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR
* OTHERWISE, ARISING FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE SOFTWARE.
*/
package ptolemaeus.math.linear.householderreductions;
import java.util.function.Consumer;
import javax.annotation.concurrent.Immutable;
import ptolemaeus.commons.Validation;
import ptolemaeus.math.commons.NumericalMethodsUtils;
import ptolemaeus.math.linear.complex.ComplexMatrix;
import ptolemaeus.math.linear.complex.ComplexVector;
import edu.umd.cs.findbugs.annotations.SuppressFBWarnings;
/** A {@link ComplexHessenbergReduction} is a reduction of a {@link ComplexMatrix matrix A} into the product
* {@code A = QRQ}H where {@link ComplexMatrix Q} is unitary and {@code R} is an {@code n}-Hessenberg matrix;
* i.e., a matrix for which all entries below the the {@code n}th sub diagonals are zero
*
*
* {@link #getR() R} is typically denoted {@code H} but, to avoid confusion with the conjugate transpose notation in
* monospaced fonts, I'm sticking with {@code R} ({@code R} due to the similarity of this to the
* {@link ComplexQRDecomposition}).
*
*
* {@code A = QRQ}H is a similarity transformation on {@code A}, meaning that {@code A} and {@code R} have
* the same Eigenvalues; i.e., Λ{@code (A) == }Λ{@code (R)}
*
*
* To compute the {@link ComplexHessenbergReduction} of a {@link ComplexMatrix} {@code A}, we proceed in much the same
* way as in a {@link ComplexQRDecomposition}, but for each iteration {@code i}, instead of computing a
* {@link Householder} matrix that "zeroes out" everything in the {@code i}-th column under the diagonal, we use
* {@link Householder} matrices to zero-out all entries in each column under the {@code n}th sub-diagonal.
*
*
* In addition, the multiplication in each iteration is different: we compute {@code A_i+1 = Q_i A_i Q_i}H,
* where the {@code Q_i} matrices are the Household rotation matrices for the ith column of {@code A_i}. The final value
* of {@code A_i} is what we save to {@code R}.
*
*
* For reference, see [1] Numerical Methods for General and Structured Eigenvalue Problems by Kressner
* (2006)
*
*
* @author Ryan Moser, [email protected] */
@Immutable
public class ComplexHessenbergReduction {
/** The Hessenberg {@link ComplexMatrix} {@code R} from {@code A = QRQ}H */
private final ComplexMatrix R;
/** the intermediate {@link ComplexMatrix} {@code Qi}s that are used to reduce the original matrix to n-Hessenberg
* form. These are cached and used to compute {@link #getQ()} as needed. Once computed, this variable is set to
* {@code null} to encourage garbage collection in case this instance is meant to be long-lived. */
private volatile ComplexMatrix[] Qs;
/** The unitary {@link ComplexMatrix} {@code Q} from {@code A = QRQ}H */
private volatile ComplexMatrix Q;
/** Constructor
*
* @param Qs the {@link ComplexMatrix matrices} {@code Q_i} whose product is equal to {@link #getQ() Q}. These are
* not saved.
* @param R {@code R} from {@code A = QRQ}H
*
* @implSpec we multiply the {@code Q} matrices together in reverse order because it results in smaller error in the
* unit tests.
*
* Intuitively, the last {@code Q} matrices in the sequence are approximately the identity matrix with a
* few elements near the lower-right that are non-zero, and including those first ensures that they "get
* the chance to contribute".
*
* @apiNote this constructor is only {@code protected} so that {@link ComplexHessenbergReduction} can be
* {@code extended} (though that seems unlikely) and for easier unit testing. */
protected ComplexHessenbergReduction(final ComplexMatrix[] Qs, final ComplexMatrix R) {
this.R = R;
this.Qs = Qs;
this.Q = null;
}
/** Constructor
*
* @param R {@code R} from {@code A = QRQ}H
* @param Q {@code Q} from {@code A = QRQ}H
*
* @apiNote this constructor is only {@code protected} so that {@link ComplexHessenbergReduction} can be
* {@code extended} (though that seems unlikely) and for easier unit testing. */
protected ComplexHessenbergReduction(final ComplexMatrix Q, final ComplexMatrix R) {
this.R = R;
this.Q = Q;
}
/** Compute a 1-{@link ComplexHessenbergReduction} of {@code A} s.t. {@code A = QRQ}H where {@code R} is
* 1-Hessenberg (all sub-diagonal entries beyond the first sub-diagonal are zero).
*
*
* E.g.:
*
*
* a a a a a a
* a a a a a a
* 0 a a a a a
* 0 0 a a a a
* 0 0 0 a a a
* 0 0 0 0 a a
*
*
* where {@code a} simply indicates a non-zero value
*
*
* From [1] Algorithm 1: Reduction to Hessenberg form (generalized w.l.o.g. to A in Cn x n):
*
*
* Input: A in Cn x n
* Output: Matrices Q and R in Cn x n where Q is Hermitian, and R is a Hessenberg matrix similar to A
*
*
* Q <- In
* R <- A
* for j <- 1, ..., n - 2, do
* H <- Hj+1(A . ej)
* Q <- Q . H
* R <- H . R . H (we can neglect the conjugate transpose because H is Hermitian)
* end for
*
*
* Where A . ej extracts the jth column of A, and Hj+1 is a function that computes
* the Householder transformation to "zero-out" the rows of A . ej after row j+1
*
*
*
*
* @param A the matrix to reduce
* @return the {@link ComplexHessenbergReduction} */
public static ComplexHessenbergReduction of(final ComplexMatrix A) {
return of(A, 1, null);
}
/** Compute an {@code n}-{@link ComplexHessenbergReduction} of {@code A} s.t. {@code A = QRQ}H where
* {@code R} is {@code n}-Hessenberg (all sub-diagonal entries beyond the first {@code n}-sub-diagonal are
* zero).
*
*
* E.g., a 2-Hessenberg matrix:
*
*
* a a a a a a
* a a a a a a
* a a a a a a
* 0 a a a a a
* 0 0 a a a a
* 0 0 0 a a a
*
*
* where {@code a} simply indicates a non-zero value
*
* @param A the matrix to reduce. Must be square (unlike in {@link ComplexQRDecomposition}).
* @param hessenbergOffset the number of sub-diagonals that we want to preserve
* @return the {@code n}-{@link ComplexHessenbergReduction}
*
* @implSpec notice how similar this implementation is to {@link ComplexQRDecomposition#of(ComplexMatrix)}. I
* considered consolidating them to reuse the relevant pieces, but I got about half-way there and things
* were becoming rather unclear, so I bit the bullet and kept them separate */
public static ComplexHessenbergReduction of(final ComplexMatrix A, final int hessenbergOffset) {
return of(A, hessenbergOffset, null);
}
/** Compute an {@code n}-{@link ComplexHessenbergReduction} of {@code A} s.t. {@code A = QRQ}H where
* {@code R} is {@code n}-Hessenberg (all sub-diagonal entries beyond the first {@code n}-sub-diagonal are
* zero).
*
*
* E.g., a 2-Hessenberg matrix:
*
*
* a a a a a a
* a a a a a a
* a a a a a a
* 0 a a a a a
* 0 0 a a a a
* 0 0 0 a a a
*
*
* where {@code a} simply indicates a non-zero value
*
* @param A the matrix to reduce. Must be square (unlike in {@link ComplexQRDecomposition}).
* @param hessenbergOffset the number of sub-diagonals that we want to preserve
* @param evaluationMonitor a {@link Consumer} the monitors the reduction. It is updated with value of
* {@link #getR() R} on each iteration.
* @return the {@code n}-{@link ComplexHessenbergReduction}
*
* @implSpec notice how similar this implementation is to {@link ComplexQRDecomposition#of(ComplexMatrix)}. I
* considered consolidating them to reuse the relevant pieces, but I got about half-way there and things
* were becoming rather unclear, so I bit the bullet and kept them separate */
public static ComplexHessenbergReduction of(final ComplexMatrix A,
final int hessenbergOffset,
final Consumer evaluationMonitor) {
return doInPlace(null, A.copy(), hessenbergOffset, evaluationMonitor);
}
/** Reduce {@code A} to {@code n}-{@link ComplexHessenbergReduction} form s.t. the original
* {@code A = QRQ}H where {@code R} is {@code n}-Hessenberg (all sub-diagonal entries beyond the first
* {@code n}-sub-diagonal are zero).
*
*
* E.g., a 2-Hessenberg matrix:
*
*
* a a a a a a
* a a a a a a
* a a a a a a
* 0 a a a a a
* 0 0 a a a a
* 0 0 0 a a a
*
*
* where {@code a} simply indicates a non-zero value
*
* @param Q the {@link ComplexMatrix} in which we will accumulate {@link #getQ()} This Matrix is not copied and
* will modified! See the {@code implNote} below. This option is useful in performance sensitive contexts
* where - e.g. - we don't want to copy the input matrices over and over again.
* @param R the matrix to reduce. Must be square (unlike in {@link ComplexQRDecomposition}). This Matrix is not
* copied and will modified! This option is useful in performance sensitive contexts where - e.g. - we
* don't want to copy the input matrices over and over again.
* @param hessenbergOffset the number of sub-diagonals that we want to preserve
* @param evaluationMonitor a {@link Consumer} the monitors the reduction. It is updated with value of
* {@link #getR() R} on each iteration.
* @return the {@code n}-{@link ComplexHessenbergReduction}
*
* @implNote notice that, when accumulating Q in-place here, that we're not accumulating Q in the same way that Q is
* lazily initialized in ::getQ when we're not updating Q in-place. When computing Q directly, we multiply
* then in the order opposite of how they're computed. Here, we multiply them in order but, for this to
* work, we have to post-multiply in-place rather than pre-multiply in-place. The practical difference is
* negligible, as documented in the tests.
*
* @implSpec notice how similar this implementation is to {@link ComplexQRDecomposition#of(ComplexMatrix)}. I
* considered consolidating them to reuse the relevant pieces, but I got about half-way there and things
* were becoming rather unclear, so I bit the bullet and kept them separate */
@SuppressWarnings("null") // using updateQInPlace prevents any NPEs, but Eclipse doesn't understand that
public static ComplexHessenbergReduction doInPlace(final ComplexMatrix Q,
final ComplexMatrix R,
final int hessenbergOffset,
final Consumer evaluationMonitor) {
NumericalMethodsUtils.requireSquare(R);
final boolean updateQInPlace = Q != null;
if (updateQInPlace) {
NumericalMethodsUtils.requireSquare(Q);
Validation.requireEqualsTo(R.getColumnDimension(), Q.getColumnDimension(), "R and Q dimensions");
}
final int dim = R.getColumnDimension();
final int n = dim - 1 - hessenbergOffset;
final ComplexMatrix[] Qs = updateQInPlace ? null : new ComplexMatrix[n];
for (int i = 0; i < n; i++) {
final ComplexVector ithColumn = new ComplexVector(R.getColumn(i)).getSubVector(i, dim - i);
final ComplexMatrix QiPrime = Householder.computeTransform(ithColumn, hessenbergOffset);
if (updateQInPlace) { // see JavaDoc
ComplexMatrix.postMultiplyInPlace(Q, QiPrime, i, i); // modifies Q data!!
}
else {
Qs[i] = QiPrime;
}
ComplexMatrix.preMultiplyInPlace(QiPrime, R, i, i); // modifies R data!!
ComplexMatrix.postMultiplyInPlace(R, QiPrime, i, i); // modifies R data!!
if (evaluationMonitor != null) {
evaluationMonitor.accept(R);
}
}
if (updateQInPlace) {
return new ComplexHessenbergReduction(Q, R);
}
else {
return new ComplexHessenbergReduction(Qs, R);
}
}
/** @return {@code Q} from the decomposition {@code A = QRQ}H. {@code Q} is unitary s.t.
* {@code Q}-1 == {@code Q}H or, equivalently, {@code Q}{@code Q}H ==
* {@code Q}H{@code Q} == {@code I}. Implicitly, {@code Q} is square. */
@SuppressFBWarnings({ "EI" })
public ComplexMatrix getQ() {
if (this.Q == null) {
synchronized (this) {
if (this.Q == null) {
final ComplexMatrix tempQ = ComplexMatrix.identity(this.R.getColumnDimension());
for (int i = this.Qs.length - 1; i >= 0; i--) {
final ComplexMatrix Qi = this.Qs[i];
ComplexMatrix.preMultiplyInPlace(Qi, tempQ, i, i); // updates tempQ!
}
this.Q = tempQ;
this.Qs = null;
}
}
}
return this.Q;
}
/** @return {@code R} from the decomposition {@code A = QRQ}H. {@code R} is {@code n}-Hessenberg, meaning
* all sub-diagonal entries beyond the first {@code n}-sub-diagonal are zero. */
@SuppressFBWarnings({ "EI" })
public ComplexMatrix getR() {
return this.R;
}
@Override
public String toString() {
return "Q:%n%s%n%nR:%n%s%n%nQ^H:%n%s".formatted(this.getQ(), this.getR(), this.getQ().conjugateTranspose());
}
}
© 2015 - 2025 Weber Informatics LLC | Privacy Policy