Estimación (lineal) cuando la covarianza del ruido es deficiente (singular)

Question

Estoy atascado con un problema que esperaba que fuera sencillo, pero estoy teniendo problemas para encontrar la manera de resolverlo con rigor. En su caso más sencillo, la situación es una estimación lineal a partir de mediciones lineales ruidosas, es decir, observamos

$$\mathbf y = \mathbf H \mathbf \theta + \mathbf w,$$

donde $\mathbf H$ es $N \times M$ para $M<N$ de rango completo (es decir, estamos en el caso sobredeterminado) y conocido, y $\mathbf w \sim {\mathcal N}(\mathbf 0, \mathbf R_w)$ es gaussiano de media cero. La cuestión es que ${\rm rank}(\mathbf R_w)=r<N$ y tengo que averiguar cómo lidiar con esto.

El enfoque estándar para este tipo de problemas sería mirar la log-verosimilitud, que da algo así como ${\rm const} - \frac 12 (\mathbf y - \mathbf H \mathbf \theta)^{\rm T}\mathbf R_w^{-1}(\mathbf y - \mathbf H \mathbf \theta)$ . La maximización de esta sobre $\theta$ da un estimador de máxima verosimilitud de la forma $$ \hat{\mathbf \theta} = \left(\mathbf H^{\rm T} \mathbf R_w^{-1} \mathbf H\right)^{-1} \mathbf H^{\rm T} \mathbf R_w^{-1} \mathbf y,$$ que coincide con el estimador de mínimos cuadrados (ponderados), también es insesgado y todo tipo de cosas bonitas.

La cuestión aquí es clara: desde que mi $\mathbf R_w$ no es invertible, nada de esto funciona. Ni siquiera puedo expresar bien la probabilidad. La densidad de $\mathbf y$ es singular, una expresión adecuada para ella necesitaría distribuciones. Sin embargo, esto no me ayudará a llegar a mi estimador de máxima verosimilitud.

Mi intuición me dice que debería dejar de proyectar $\mathbf y$ en un espacio de menor dimensión donde el vector de ruido efectivo tiene entonces rango completo. Pero no estoy seguro de cómo hacer esto y demostrar que la estimación resultante sigue siendo óptima. ¿O la respuesta es que un estimador ML no existe en este caso y por lo tanto no está claro lo que significa la optimalidad? Eso me sorprendería. Después de todo, este tipo de cosas pueden ocurrir fácilmente, basta con repetir una de mis observaciones en $\mathbf y$ . Espero una respuesta sencilla, que no conduzca a una mejor estimación.

Estoy pensando que el camino a seguir es definir alguna dimensión inferior $\mathbf z = \mathbf P \mathbf y = \mathbf P \mathbf H \mathbf \theta + \mathbf P \mathbf w$ donde el vector de ruido efectivo $\mathbf P \mathbf w$ tiene una covarianza de rango completo y luego demostrar que $\mathbf z$ es una estadística suficiente para estimar $\theta$ . Pero, ¿cómo construyo $\mathbf P$ ? Empecé a tomar el $r$ vector propio dominante de $\mathbf R_w$ pero esto me parece mal: si me proyecto en este espacio, estoy matando una parte de $\mathbf H$ puede ocurrir que $\mathbf H$ vive en el subespacio ortogonal en parte, o incluso en su totalidad. Esto último es curioso ya que en este caso podría proyectar $\mathbf y$ en el espacio de la columna de $\mathbf H$ y obtendría una observación libre de ruido tal que pueda estimar $\mathbf \theta$ exactamente. Esto podría significar que necesito algunas condiciones en el rango de $\mathbf H$ y $\mathbf R_w$ .

Aun así, intuitivamente debería haber una solución sencilla que coincida con el WLS para el caso especial $r=N$ . ¿Puede alguien ayudarme a arrojar algo de luz sobre esto?

Answer 1

Después de haber reflexionado un poco más sobre esto, ahora puedo responder a mi propia pregunta. Lo más importante y posiblemente útil para otros: en el caso $\mathbf R_{\mathbf w}\succeq \mathbf 0$ la probabilidad puede escribirse como

$$ f_{\mathbf y}(\mathbf y | \mathbf \theta) = \begin{cases} {\rm const} \cdot {\rm e}^{-\frac 12 (\mathbf y - \mathbf H \mathbf \theta)^{\rm T} \mathbf R_{\mathbf w}^+ (\mathbf y - \mathbf H \mathbf \theta)} & \mathbf R_{\mathbf w}^\perp \cdot(\mathbf y - \mathbf H \mathbf \theta) = \mathbf 0 \\ 0 & \text{otherwise}, \end{cases}$$

donde $\mathbf R_{\mathbf w}^+$ es una inversa generalizada y $\mathbf R_{\mathbf w}^\perp$ es un proyector sobre el complemento de $\mathbf R_{\mathbf w}$ . La mejor manera de tratar la restricción es descomponer $\mathbf \theta$ en un vector en el rango de $\mathbf R_{\mathbf w}^\perp \mathbf H$ que satisface la restricción y otro vector en su complemento que puede ser elegido libremente. A continuación, podemos formar la derivada de la (log-)verosimilitud con respecto a este último vector y ponerla a cero para encontrar el estimador ML.

Como se esperaba, el estimador ML resultante es lineal, insesgado y alcanza el límite de Cramer-Rao. Por tanto, también es AZUL. Como ventaja adicional, para el caso en que $r < N-M$ podemos satisfacer las restricciones en el espacio del complemento de $\mathbf R_{\mathbf w}$ y conseguir un estimador con varianza cero, como se esperaba.

editar Ya que me pidieron el estimador ML: No he encontrado una forma súper corta de escribirlo, aquí está mi mejor intento actual: $$ \begin{align} \mathbf{\hat{\theta}}_{\rm ML} = \left(\mathbf F^+ \mathbf P^{\rm T} + \mathbf F^\perp \left[\mathbf F^{\perp^{\rm T}} \mathbf H^{\rm T} \mathbf \Sigma^+ \mathbf H \mathbf F^\perp\right]^{-1} \mathbf F^{\perp^{\rm T}} \mathbf H^{\rm T} \mathbf \Sigma^+ \left(\mathbf I - \mathbf H \mathbf F^+ \mathbf P^{\rm T} \right)\right)\mathbf y. \end{align} $$ donde $\mathbf F = \mathbf P^{\rm T} \mathbf H\in \mathbb{R}^{N-r \times M}$ y $\mathbf P \in \mathbb{R}^{N \times N-r}$ es una base ortonormal para el espacio nulo de $\mathbf R_w$ .

Tenga en cuenta que para $r < N-M$ este estimador se simplifica a $$\mathbf{\hat{\theta}}_{\rm ML} = \mathbf F^+ \mathbf P^{\rm T} \mathbf y.$$

En este caso, la inserción de $\mathbf y = \mathbf H \mathbf \theta + \mathbf w$ da $\mathbf{\hat{\theta}}_{\rm ML} = \mathbf \theta$ (ya que $\mathbf P^{\rm T} \mathbf w = \mathbf 0$ ), es decir, obtenemos una estimación exacta con varianza cero.