Why isn't the pseudo-inverse the best choice in my linear estimation problem?

Question

Context:

I have a problem of the following form. Let $\boldsymbol{\theta}\in\mathbb{R}^n$ be a fixed vector I want to estimate. Let $\mathbf{M}\in\mathbb{R}^{m\times n}$ be a matrix with $m>n$ and $\text{rank}(\mathbf{M})=n$. Moreover, let $\mathbf{w}\in\mathbb{R}^m$ be a vector with independent identically distributed components, each distributed as $\sim\mathcal{N}(0,\sigma^2)$. Assume that you have access to the "noisy vector": $$ \mathbf{x} = \mathbf{M}\boldsymbol{\theta} + \mathbf{w} $$ The question is to obtain the "best" estimate $\hat{\boldsymbol{\theta}}$ of $\boldsymbol{\theta}$, under mean squared error criterion. My first instinct is to set my estimate using a pseudo-inverse according to $$ \hat{\boldsymbol{\theta}}_1 = (\mathbf{M}^\top\mathbf{M} )^{-1}\mathbf{M}^\top\mathbf{x} $$

However, in a particular problem I have, I found by "accident" another matrix $\mathbf{M}_{\text{special}}$ which seems to perform better than $(\mathbf{M}^\top\mathbf{M} )^{-1}\mathbf{M}^\top$ somehow. Call $$ \hat{\boldsymbol{\theta}}_2 = \mathbf{M}_{\text{special}}\mathbf{x} $$ for this special matrix I found. The concrete data I'm using can be found by running the matlab script at the end of this question.

To test which estimate performs better, I'm computing the "overall" errors $\sqrt{\frac{1}{n}({\boldsymbol{\theta}}-\hat{\boldsymbol{\theta}}_1)^\top({\boldsymbol{\theta}}-\hat{\boldsymbol{\theta}}_1)}$ and $\sqrt{\frac{1}{n}({\boldsymbol{\theta}}-\hat{\boldsymbol{\theta}}_2)^\top({\boldsymbol{\theta}}-\hat{\boldsymbol{\theta}}_2)}$ and comparing them. My observations are the following:

• If $\sigma^2=0$ (no noise), then I can recover $\hat{\boldsymbol{\theta}}_1={\boldsymbol{\theta}}$, which is expected since $\mathbf{M}^\top\mathbf{M}$ is full rank. On the other hand ${\boldsymbol{\theta}}-\hat{\boldsymbol{\theta}}_2$ is a non-zero vector (even without noise). Hence, the pseudo-inverse estimate performs better.
• As I increase the noise variance, both estimates decrease performance (which is expected). However, the pseudo-inverse one decreases faster. In fact, when $\sigma=0.5$, both estimates perform roughly the same.
• For bigger values of $\sigma^2>1$, I find that $\hat{\boldsymbol{\theta}}_2$ ALWAYS performs better than $\hat{\boldsymbol{\theta}}_1$.

Question:

How is it possible for a different linear estimate $\hat{\boldsymbol{\theta}}_2$ outperform what is supposed to be the best linear estimator $\hat{\boldsymbol{\theta}}_1$? I assume I misunderstood something about the "optimality nature" of the pseudo-inverse estimate I'm using. However, I cannot find the error in my reasoning. What makes $\hat{\boldsymbol{\theta}}_2$ better?

Data:

The data and script I'm using can be found here: https://github.com/RodrigoAldana/se_question_data

I know that I'm generating new $\boldsymbol{\theta}$ in each run. Please note that I also tried to fix $\boldsymbol{\theta}$ for all runs obtaining similar results.

Answer 1

$(M'M)^{-1}M'x$ is the best linear unbiased estimator and has a covariance matrix $\sigma^2(M′M)^{−1}$, which differs from all other unbiased linear estimators by a positive semidefinite matrix. Your estimator is biased, possibly among other things, so the Gauss-Markov theorem does not apply to it.

It is well-known that biased linear estimators, e.g., Ridge regression estimators, can outperform OLS. Perhaps you found one such!