# 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$$.

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}$$