How to compute an indefinite generalisation of QR decomposition

Question

Given an arbitrary complex matrix $M$ and real, diagonal but possibly indefinite matrix $\Delta$, the problem is to solve the following system of equations: $$\begin{aligned} M^*\Delta M &= LD^2L^*\\ M &= QDL^* \end{aligned}$$ for lower unitriangular $L$, diagonal $D$, and $Q$. The running time should be measured in terms of the number of multiplications.

Notice that when $\Delta$ is the identity matrix, this problem is exactly equivalent to QDR decomposition, and can be solved using any algorithm for the latter. The cost is $4n^3/3$. (If you're unfamiliar with QDR decomposition, it's essentially QR decomposition followed by factorising $R$ into $DR'$ where $D$ is diagonal and $R'$ is unitriangular).

When $\Delta$ is a positive diagonal matrix, the problem can again be reduced to QDR decomposition. The cost is therefore again $4n^3/3$.

When $\Delta$ is indefinite, it's not clear if a $4n^3/3$ algorithm exists. One can proceed by multiplying out the matrices $M^* \Delta M$ (which costs $n^3$), taking the LDL decomposition using the Cholesky-Crout algorithm (which costs $n^3/3$) and then finding $Q$ by forward substitution (which costs $n^3/2$). The overall cost is thus $11n^3/6 > 4n^3/3$.

I'm wondering if there's a way to generalise Gram-Schmidt, or any other algorithm for QDR decomposition, so that in the indefinite case the cost is still $4n^3/3$.

Answer 1

The problem can be solved by recognising that it's essentially the same as QR decomposition, but where the role of the inner product is replaced by the diagonal sesquilinear form $\langle u, v \rangle = \sum_{i=1}^n \delta_i \overline u_i v_i$. Therefore one can take the Gram-Schmidt algorithm from Wikipedia and change the role of the inner product to the above sesquilinear form.

One must make one further modification: The expression for $R$ given on Wikipedia is not valid for an arbitrary diagonal sesquilinear form: $$R = \begin{pmatrix} \langle\mathbf{e}_1, \mathbf{a}_1\rangle & \langle\mathbf{e}_1, \mathbf{a}_2\rangle & \langle\mathbf{e}_1, \mathbf{a}_3\rangle & \ldots \\ 0 & \langle\mathbf{e}_2, \mathbf{a}_2\rangle & \langle\mathbf{e}_2, \mathbf{a}_3\rangle & \ldots \\ 0 & 0 & \langle\mathbf{e}_3, \mathbf{a}_3\rangle & \ldots \\ \vdots & \vdots & \vdots & \ddots \end{pmatrix}$$ Rather, it must be $$R = \begin{pmatrix} \frac{\delta_1}{|\delta_1|}\langle\mathbf{e}_1, \mathbf{a}_1\rangle & \frac{\delta_1}{|\delta_1|}\langle\mathbf{e}_1, \mathbf{a}_2\rangle & \frac{\delta_1}{|\delta_1|}\langle\mathbf{e}_1, \mathbf{a}_3\rangle & \ldots \\ 0 & \frac{\delta_2}{|\delta_2|}\langle\mathbf{e}_2, \mathbf{a}_2\rangle & \frac{\delta_2}{|\delta_2|}\langle\mathbf{e}_2, \mathbf{a}_3\rangle & \ldots \\ 0 & 0 & \frac{\delta_3}{|\delta_3|}\langle\mathbf{e}_3, \mathbf{a}_3\rangle & \ldots \\ \vdots & \vdots & \vdots & \ddots \end{pmatrix}$$