# The question

Is there a variant of the Householder QR algorithm to orthonormalize a set of vectors with respect to an inner product if no orthonormal basis is known a priori?

# Background

Let's assume that $A\in\mathbb{C}^{n,k}$ is a matrix with $\operatorname{rank}(A)=k$. We can then compute the QR factorization $A=QR$ with $Q=[q_1,\ldots,q_k]\in\mathbb{C}^{n,k}$ and $R\in\mathbb{C}^{k,k}$ to obtain an orthonormal basis (with respect to the Euclidean inner product) $q_1,\ldots,q_k$ of $\operatorname{range}(A)$, i.e. $q_i^*q_j = \delta_{ij}$.

If orthonormality is critical in the presence of round-off errors, the Householder QR algorithm is the method of choice (cf. Golub, Van Loan. Matrix Computations.). Compared to the (modified) Gram-Schmidt algorithm the Householder QR algorithm benefits from the favorable round-off properties of Householder transformations.

I am aware of the paper Trefethen, Householder triangularization of a quasimatrix, 2009 where the $L^2$-inner product is used for a Householder QR algorithm. However, an orthonormal basis (the Legendre polynomials) is known in the paper by Trefethen.

In my application only the Hermitian and positive definite matrix $B\in\mathbb{C}^{n,n}$ defining the inner product $\langle x,y \rangle_B = x^* B y$ is known. Because $n\gg k$ it's also infeasible to compute an eigen-decomposition of B.

It seems to me that knowing one orthonormal basis (i.e. the standard basis) is crucial in the construction of the Householder QR algorithm but perhaps I'm wrong and there is a trick to let Householder QR work with arbitrary inner products?

I haven't read Trefethen's paper and I don't have access to journals at the moment. However, I think you might find these two papers interesting:

1. Singer, "Indefinite QR factorization"
2. Singer, Singer, "Orthosymmetric block reflectors",

Some other papers, either theirs or just Sanja Singer's, might be relevant as well.

These revolve around the so called orthosymmetric scalar products, induced by a matrix $J$ (instead of your $B$). "Orthosymmetric" means that $J^* = \tau J$, where $|\tau| = 1$. Your $B$ is a very, very special case of such matrix (usually, Hermitian $J$ is the most interesting special case).

The worse obstacle in the above papers is the one you don't have: generally, $x^* J x$ can be zero or even negative, but your $B$ is positive definite, so this will never happen to you. This means that you can skip any mentions of the "degenerate vectors" in the above papers, while the prinicples and ideas should still work.

Most of the [1] revolves around the scalar products induced by a signature matrix $J = \mathop{\rm diag}(j_1,\dots,j_n)$, where $j_k \in \{-1,1\}$. This means that $J$ is diagonal, which your $B$ is not. However, if I remember that paper right, there should be some mention of the more general case where $J$ is any Hermitian matrix. Also, in [2], the authors deal with the (block) reflectors in general orthosymmetric scalar product spaces, which is far more general than you need, so it should help you solve your problem.

• Thanks for your answer and the references. The construction of Householder transformations in a non-Euclidean inner product is not a problem. The actual problem is that you don't have a basis at hand that is orthonormal with respect to the given inner product. So you don't know the vectors your Householder transformation should map to. The question still is if it's possible to give an applicable algorithm. Unfortunately, the paper by Singer and Singer with the promising title Householder-like algorithms for orthosymmetric QR factorizations (cited in [1]) apparently was never published. :( – andrenarchy Jul 2 '13 at 9:39
• You are right, it is still a work in progress, but you might be O.K. with chapter 4.3. in Sanja's doctoral dissertation. Unfortunately, it is in Croatian, so you'll need to find someone to translate it for you, or try to deduce the ideas from the formulas. – Vedran Šego Jul 3 '13 at 13:32
• Unfortunately I don't understand Croatian and I don't have anybody around who could translate. Automatic translation fails completely (as expected) ;) so I'd still leave this question open, perhaps there are other ideas (or translations) out there! :) – andrenarchy Jul 3 '13 at 21:33