QL decomposition

 A couple days ago I encountered a really neat solution to creating a basis using a Cholesky decomposition. Particularly, when we have a set of basis function $Y = {y_0,y_1,\cdots , y_n}$, we will want a set of coefficient vectors  $C = {\vec{c}_0,\vec{c}_1,\cdots , \vec{c}_n}$ such that $\langle \sum_j^n c_{k,j} y_j,  \sum_i^n c_{l,i} y_i\rangle = \delta_{k,l}$, and, ideally we also wish the constraint $c_{k,j} = 0 \forall j>k$. More succinctly, an orthonormal basis where there are only order $k$ or lower terms for the $k$-th basis element. Putting that into math terms (and using some abusive notation for the matrix of inner products) we want to find $R$ such that $R^TY^TYR = I$.

The thinking man would use a cholesky decomposition and an inverse to say that if $Y^TY = LL^T$ then $L^{-1}Y^TYL^{-T} = I$, and as $L^{-T}$ is upper triangular, we're finished.

I however am not a thinking man, just a mindless automaton, driven by a basic response to stimuli and only distantly related to chordate as a whole. As such I only have two tools at my disposal, an SVD, and a QR decomposition. How then may I solve this?

Given symmetry we know that if $Y^TY = U \Sigma V^T$ then $U = V$ and thus $Y^TY = U \Sigma U^T$. Multiplying out on both sides we get $\Sigma^{-1/2}U^TY^TYU \Sigma^{-1/2}= I$. This is nice but not upper triangular. We note here that if we found $Q$ and $R$ such that $U \Sigma^{-1/2}=RQ$, we could then multiply $Q^TR^TY^TYRQ$ by $Q$ on the left and $Q^T$ on the right to get $R^TY^TYR = QQ^TR^TY^TYRQQ^T = QQ^T = I$.

How then do we find $R$ and $Q$? Well, we see this is equivalent to finding $Q$ and $R$ such that $\Sigma^{-1/2}U^T = Q^T R^T$. This is a QR decomposition with a lower and not upper triangular matrix.

Because common sense can not and will not prevail we do the obvious thing, which is look at what happens when we 'flip' a matrix. Formally, $flip(A^{n \times m},0)_{i,j} = A_{n-i,j}$ and $flip(A^{n \times m},1)_{i,j} = A_{i,m-j}$, that is, flipping the matrix left to right and up and down.

When we apply this flip operation to a matrix product, we see $flip(AB,0) = flip(A,0) B$ and $flip(AB,1) = A flip(B,1)$.

Suppose then that we took a QR of $flip(\Sigma^{-1/2}U^T,1) = \tilde{Q} \tilde{R}$, we see then that $\tilde{Q} ^T flip(\Sigma^{-1/2}U^T,1) = \tilde{R}$, and then that $flip(\tilde{Q} ^T flip(\Sigma^{-1/2}U^T,1),0) = flip(\tilde{Q} ^T,0) flip(\Sigma^{-1/2}U^T,1) = flip(\tilde{R},0)$, and then $flip(flip(\tilde{Q} ^T,0) flip(\Sigma^{-1/2}U^T,1),1) =flip(\tilde{Q} ^T,0) \Sigma^{-1/2}U^T = flip(flip(\tilde{R},0),1)$. Thus, as $flip(flip(\tilde{R},0),1)$ is lower triangular we get a solution for $\Sigma^{-1/2}U^T = Q^T R^T$ where $R = flip(flip(\tilde{R},0),1)^T$. Here's some code you can run to verify it for yourself.

u,s,v = np.linalg.svd(Y.T @ Y)

c = (u/np.sqrt(s))

q,r = np.linalg.qr(np.flip(c.T,axis=1));

c = np.flip(np.flip(r,axis=0),axis=1).T

This process can be used to generate decompositions for every form and order of triangular matrix by orthonormal matrix.