How to compute an indefinite generalisation of QR decomposition 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$.
 A: 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}$$
