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


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


Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.