[proofplan] We compare the variance of the scalar linear estimator $c^\top AY$ with that of the ordinary least squares estimator $c^\top BY$. Since both $A$ and $B$ are unbiased estimator matrices, the matrix $A-B$ annihilates the column space of $X$, which makes the cross term in the variance expansion vanish. The variance difference is therefore exactly $\sigma^2 |(A-B)^\top c|^2$, and the equality condition follows from positivity of $\sigma^2$ and the Euclidean norm. [/proofplan] [step:Verify that both estimator matrices reproduce $\beta$ on the model space] Since $X$ has full column rank, $X^\top X \in \mathbb{R}^{p \times p}$ is invertible. Therefore $B=(X^\top X)^{-1}X^\top$ is well-defined, and \begin{align*} BX = (X^\top X)^{-1}X^\top X = I_p. \end{align*} By hypothesis, $AX=I_p$. Hence \begin{align*} (A-B)X = AX-BX = I_p-I_p=0. \end{align*} Equivalently, \begin{align*} X^\top(A-B)^\top=0. \end{align*} [/step] [step:Expand the variance difference and cancel the cross term] Fix $c \in \mathbb{R}^p$. Define the vectors \begin{align*} u &:= B^\top c \in \mathbb{R}^n, & v &:= (A-B)^\top c \in \mathbb{R}^n. \end{align*} Then \begin{align*} A^\top c = B^\top c + (A-B)^\top c = u+v. \end{align*} For any fixed vector $w \in \mathbb{R}^n$, the scalar random variable $w^\top Y$ has variance \begin{align*} \operatorname{Var}(w^\top Y) &= \operatorname{Var}(w^\top(X\beta+\varepsilon)) \\ &= \operatorname{Var}(w^\top \varepsilon) \\ &= w^\top \operatorname{Var}(\varepsilon) w \\ &= \sigma^2 w^\top w = \sigma^2 |w|^2. \end{align*} Applying this with $w=A^\top c$ and $w=B^\top c$ gives \begin{align*} \operatorname{Var}(c^\top AY)-\operatorname{Var}(c^\top BY) &= \sigma^2 |A^\top c|^2-\sigma^2 |B^\top c|^2 \\ &= \sigma^2 |u+v|^2-\sigma^2 |u|^2 \\ &= \sigma^2\left(2u^\top v+|v|^2\right). \end{align*} It remains to show $u^\top v=0$. Using the definitions of $u$ and $v$, \begin{align*} u^\top v &= (B^\top c)^\top (A-B)^\top c \\ &= c^\top B(A-B)^\top c \\ &= c^\top (X^\top X)^{-1}X^\top(A-B)^\top c \\ &= 0, \end{align*} because $X^\top(A-B)^\top=0$. Therefore \begin{align*} \operatorname{Var}(c^\top AY)-\operatorname{Var}(c^\top BY) = \sigma^2 |(A-B)^\top c|^2. \end{align*} Since $\hat{\beta}=BY$, this is the same as \begin{align*} \operatorname{Var}(c^\top AY)-\operatorname{Var}(c^\top \hat{\beta}) = \sigma^2 |(A-B)^\top c|^2. \end{align*} [guided] Fix $c \in \mathbb{R}^p$. The scalar estimators $c^\top AY$ and $c^\top BY$ are obtained by applying the row vectors $c^\top A$ and $c^\top B$ to $Y$. It is therefore convenient to transpose these row vectors and work with vectors in $\mathbb{R}^n$. Define \begin{align*} u &:= B^\top c \in \mathbb{R}^n, & v &:= (A-B)^\top c \in \mathbb{R}^n. \end{align*} Then $A^\top c=u+v$. We first compute the variance of an arbitrary scalar linear form in $Y$. Let $w \in \mathbb{R}^n$. Since $Y=X\beta+\varepsilon$ and $w^\top X\beta$ is deterministic, \begin{align*} \operatorname{Var}(w^\top Y) &= \operatorname{Var}(w^\top X\beta+w^\top\varepsilon) \\ &= \operatorname{Var}(w^\top\varepsilon). \end{align*} The covariance matrix of $\varepsilon$ is $\sigma^2 I_n$, so the variance of the scalar random variable $w^\top\varepsilon$ is \begin{align*} \operatorname{Var}(w^\top\varepsilon) = w^\top \operatorname{Var}(\varepsilon)w = w^\top(\sigma^2 I_n)w = \sigma^2 |w|^2. \end{align*} Applying this formula first to $w=A^\top c$ and then to $w=B^\top c$ gives \begin{align*} \operatorname{Var}(c^\top AY)-\operatorname{Var}(c^\top BY) &= \sigma^2 |A^\top c|^2-\sigma^2 |B^\top c|^2 \\ &= \sigma^2 |u+v|^2-\sigma^2 |u|^2 \\ &= \sigma^2\left(2u^\top v+|v|^2\right). \end{align*} The key point is that the cross term $u^\top v$ vanishes. Indeed, \begin{align*} u^\top v &= (B^\top c)^\top (A-B)^\top c \\ &= c^\top B(A-B)^\top c \\ &= c^\top (X^\top X)^{-1}X^\top(A-B)^\top c. \end{align*} From $AX=I_p$ and $BX=I_p$, we have $(A-B)X=0$, and transposing gives $X^\top(A-B)^\top=0$. Hence \begin{align*} u^\top v=0. \end{align*} Substituting this into the variance expansion yields \begin{align*} \operatorname{Var}(c^\top AY)-\operatorname{Var}(c^\top BY) = \sigma^2 |(A-B)^\top c|^2. \end{align*} Since $\hat{\beta}=BY$, this is exactly \begin{align*} \operatorname{Var}(c^\top AY)-\operatorname{Var}(c^\top \hat{\beta}) = \sigma^2 |(A-B)^\top c|^2. \end{align*} [/guided] [/step] [step:Read off the equality condition from the squared norm] Because $\sigma^2>0$ and $|z|^2=0$ for a vector $z \in \mathbb{R}^n$ if and only if $z=0$, the identity \begin{align*} \operatorname{Var}(c^\top AY)-\operatorname{Var}(c^\top \hat{\beta}) = \sigma^2 |(A-B)^\top c|^2 \end{align*} implies \begin{align*} \operatorname{Var}(c^\top AY)=\operatorname{Var}(c^\top \hat{\beta}) \end{align*} if and only if \begin{align*} (A-B)^\top c=0. \end{align*} [/step] [step:Use equality for all directions to identify the estimator matrix] Assume now that \begin{align*} (A-B)^\top c=0 \end{align*} for every $c \in \mathbb{R}^p$. Let $D:=A-B \in \mathbb{R}^{p \times n}$. For each column index $j \in \{1,\dots,n\}$, let $d_j \in \mathbb{R}^p$ denote the $j$-th column of $D$. The $j$-th component of $D^\top c$ is $d_j^\top c$, so the hypothesis gives \begin{align*} d_j^\top c=0 \end{align*} for every $c \in \mathbb{R}^p$. Taking $c=d_j$ gives \begin{align*} |d_j|^2=0, \end{align*} hence $d_j=0$. Since this holds for every $j \in \{1,\dots,n\}$, every column of $D$ is zero. Thus $D=0$, so $A-B=0$ and therefore $A=B$. [/step]