1

Suppose you observe some matrix $\Sigma$ and you know it is of the form $$\Sigma = H\sigma \sigma^T H^T,$$ where $H^TH=I$ and $HH^T=P$ satisfies $P^2=P$, $P^T=P$ (i.e. it is an orthogonal projection), $H\in \mathbb{R}^{D\times d}$, $\sigma\in \mathbb{R}^{d\times d}$, but $\sigma$ and $\sigma\sigma^T$ are unknown. We want to estimate $H$. Are there any standard methods for doing so?

If $\sigma\sigma^T=I$, we can just use SVD. Since in that case $\Sigma = P=HH^T$ if we perform SVD we get $\Sigma = USV^T$ and assuming the singular values are in increasing order, the first $K<D$ entries of the diagonal of $S$ are zero, so if $A=U\sqrt{S}$ and the columns are $\vec{a}_1,\dotsc, \vec{a}_D$, we take the matrix $$H = [\vec{a}_K, \vec{a}_{K+1}, \dotsc, \vec{a}_D].$$

I suspect in the case $\sigma\sigma^T\neq I$, this is too much to ask for without knowing/observing $\sigma$ or $\sigma\sigma^T$. Seems like the best we can do is try to recover $B = H\sigma$. If we know $\sigma \sigma^T$ and it is invertible, then we could get $B\sigma^T = H\sigma\sigma^T$ so that $H = B\sigma^T(\sigma\sigma^T)^{-1}$.

It might be a naive question hoping for too much but I figured it can't hurt to ask for advice. I apologies if this is something basic, and would love references for these sorts of things in linear algebra/matrix factorization.

Nap D. Lover
  • 1,207
  • 1
    One approach is to use a "spectral decomposition" of $P$, which should in a sense coincide with its SVD because $P$ is positive semidefinite – Ben Grossmann Nov 07 '23 at 18:54
  • 1
    So actually in the case of $\sigma\sigma^T \neq I$, as long as $\sigma$ is well-conditioned (singular values are close in magnitude), your approach mostly works; just take the first $d$ columns of $U$ without multiplying by $\sqrt{S}$. – Ben Grossmann Nov 07 '23 at 18:56
  • @BenGrossmann Ah, interesting! If you would like to expand this into a little answer, I would definitely accept it. – Nap D. Lover Nov 07 '23 at 19:03
  • If you're using Python or Matlab, I could give you a code snippet as well – Ben Grossmann Nov 07 '23 at 19:06
  • @BenGrossmann Oh wow, thanks! I am indeed using Python. – Nap D. Lover Nov 07 '23 at 19:08
  • 1
    Quick caveat: if $\sigma$ has repeated singular values, we can't guarantee that $H$ is recovered exactly, only that $HH^T = P$ is correct – Ben Grossmann Nov 07 '23 at 19:17
  • @BenGrossmann That is actually perfectly fine for my intentions which were not explicitly written in the OP, apologies--I have a model orthogonal projection $P_\theta$ that I want to get close to the true $P$ ("hidden" in the observed $\Sigma$). So, indeed I only really care that $HH^T=P$ is recovered correctly. – Nap D. Lover Nov 07 '23 at 19:20
  • 1
    Interesting fact - the Eckart-Young-Mirsky theorem ensures that, if $\Sigma$ is perturbed with noise, the first columns from the $U$ of the SVD give the "best" approximation by "reasonable metrics". – Ben Grossmann Nov 07 '23 at 19:31

1 Answers1

1

As noted in the comments, it suffices to extract the first $d$ columns of $U$ from the SVD $\Sigma = U\Sigma V^T$. Here's some Python code illustrating this approach.

import numpy as np
import numpy.linalg as la
from numpy.random import randn

parameters


d = 3
D = 10


generate example


sigma = randn(d,d)
H = la.qr(randn(D,d))[0]
Sigma = H @ sigma @ sigma.T @ H.T
P = H @ H.T


apply method


H = la.svd(Sigma)[0][:,:d]


check accuracy of corresponding P


print(f"Relative error: {la.norm([email protected] - P)/la.norm(P):.3g}")
Ben Grossmann
  • 225,327