Proof of the multivariate Cramer-Rao inequality I search a detailed proof of the multivariate Cramer-Rao inequality in the general case where the estimator is not necessarily unbiased.
Let $T(X)$ be an estimator of the parameter $\theta\in\mathbb{R}^m$.
Let $\psi$ be the expectation of $T(X)$:
\begin{align}
\psi(\theta) = \mathbb{E}(T(X))
\end{align}
for any $\theta\in\mathbb{R}^m$.
Assume that the following regularity condition holds:
\begin{align}
\frac{\partial }{\partial \theta} \int_{\mathbb{R}^n} T(x) f(x ; \theta) dx
= \int_{\mathbb{R}^n} \frac{\partial }{\partial \theta} T(x) f(x ; \theta) dx
\end{align}
for any $\theta\in\mathbb{R}^m$.
If $\psi$ is differentiable, therefore,
\begin{align}
\mathbf{Cov}(T(X)) \geq 
\left[ \frac{\partial \psi(\theta)}{\partial \theta} \right]^T
[I(\theta)]^{-1}
\left[ \frac{\partial \psi(\theta)}{\partial \theta} \right]
\end{align}
for any $\theta\in\mathbb{R}^m$.
The Cramer-Rao inequality is a matrix inequality.
If $A$ and $B$ are two real square matrices, we say that $A \geq B$
if the matrix $A - B$ is positive definite.
The matrix $\frac{\partial \psi(\theta)}{\partial \theta}$ is the
Jacobian matrix of $\psi$:
$$
\left(\frac{\partial \psi(\theta)}{\partial \theta}\right)_{ij}
 = \frac{\partial \psi(\theta)_i}{\partial \theta_j}
$$
for $i,j = 1, \ldots, m$.
Let $s$ be the score function, i.e. the gradient vector of the log-likelihood.
It is not very difficult to see that:
$$
\mathbf{Cov}(s_j(\theta),T(X)) = \frac{\partial \psi(\theta)_i}{\partial \theta_j}
$$
for $j=1,\ldots, m$. The previous equation can be obtained based on the straightforward differentiation of the log-likelihood, the regularity condition and the definition of $\psi$.
At this point, the univariate proof uses the Cauchy-Schwartz probabilistic inequality, but this cannot be used in the multivariate case. How to conclude?
Indeed, the univariate proof is not too difficult, but the multivariate proof is more tricky. In Jun Shao "Mathematical Statistics", p.169, the proof is "left to the reader"...
 A: I have found a helpful inequality by extending the Cauchy-Schwartz Inequality to the multivariate case. The inequality is presented in https://benmoran.wordpress.com/2016/07/16/cauchy-schwarz-for-outer-products-as-a-matrix-inequality. Although the question was asked a long time ago, I want to show the whole proof since I myself also struggled for a long time looking for an answer.

We first present the inequality we are going to use, which I called the multivariate Cauchy-Schwartz inequality. (I don't know if it has a special name.)
Lemma. (Multivariate Cauchy-Schwartz Inequality).
Suppose $y\in\mathbb{R}^m$ and $z\in\mathbb{R}^n$ are random vectors satisfying $\text{Var}(y)>0$, then $$\text{Var}(z)\ge \text{Cov}(z, y)\text{Var}(y)^{-1} \text{Cov}(y, z).$$
Proof.
Let $u:=y-\mathbb Ey$ and $v:=z-\mathbb Ez$ so that $\mathbb Eu=\mathbb Ev=0$.
For any matrix $A\in\mathbb{R}^{n\times m}$, we have $(v+Au)(v+Au)^T\ge 0$. It follows that $\mathbb E[(v+Au)(v+Au)^T]\ge 0$. Hence by linearity of expectation we have $$\mathbb E[vv^T]+A\mathbb E[uv^T]+\mathbb E[vu^T]A^T+A\mathbb E[uu^T]A^T\ge0,$$and thus by definition $$\text{Var}(vv^T)+A\text{Cov}(u, v)+\text{Cov}(v, u)A^T+A\text{Var}(u)A^T\ge0.$$
Now let $A=-\text{Cov}(v, u)\text{Var}(u)^{-1},$ ($\text{Var}(u)$ is invertible since it is positive semi-definite.) then we can eliminate the third and the last term, getting $$\text{Var}(v)-\text{Cov}(v,u)\text{Var}(u)^{-1} \text{Cov}(u,v)\ge 0$$ and the conclusion follows. Q.E.D.
Based on this inequality, we are now prepared to derive the multivariate version of the Cramer-Rao lower bound using a very similar approach as in the univariate case.
Proof of the Cramer-Rao lower bound.
Suppose $T$'s domain is $\mathbb R^k$. ($k=m$ if $T(X)$ is an estimator of $\theta$, and I separate them because it can help clarify the confusing matrix calculus part.) Let $y=\frac{\partial}{\partial \theta} \log f(X; \theta)\in \mathbb R^m$ and $z=T(X)\in \mathbb R^k$. We give 2 claims as follows.
Claim 1. $I(\theta)=\text{Var}(y)$.
Proof. The proof is very similar to the univariate case. Given that
$$
\begin{align*}
\mathbb Ey&:=\mathbb E[\frac{\partial}{\partial \theta} \log f(X; \theta)]\\
&=\int \frac{\partial}{\partial \theta} \log f(X; \theta) f(X; \theta)d\nu\\
&=\int \frac{\partial}{\partial \theta} f(X; \theta)d\nu\\
&=\frac{\partial}{\partial \theta}\int f(X; \theta)d\nu=0,\end{align*}$$ we have $$\text{Var}(y)=\mathbb E[(y-\mathbb Ey)(y-\mathbb Ey)^T]=\mathbb E[yy^T]=I(\theta).$$
Claim 2. $\frac{\partial}{\partial \theta} \psi(\theta)=\text{Cov}(y, z)$.
Proof. This part is also straightforward.
$$
\begin{align*}
\frac{\partial}{\partial \theta} \psi(\theta)&=\frac{\partial}{\partial \theta}\int T(X)f(X; \theta)d\nu\\
&=\int \frac{\partial}{\partial \theta} f(X; \theta) T(X)^Td\nu\\
&=\int \frac{\partial}{\partial \theta} \log f(X; \theta) T(X)^T f(X; \theta)d\nu\\
&=\mathbb E[yz^T]\\
&=\text{Cov}(y, z). 
\end{align*}$$
One thing worth noting is that here we transpose $T(X)$ when applying the regularity condition and putting the partial derivative inside, in order for the matrix shape to agree with each other. It is a bit different from the Jacobian, since $I(\theta)\in\mathbb R^{m\times m}$ and we want $\frac{\partial}{\partial \theta} \psi(\theta)$ to be of shape $m\times k$. (If there is anything wrong, please correct me.)
Based on the 2 claims and the inequality, the statement we want to prove emerges naturally:
$$\text{Var}(T(x))=\text{Var}(z)\ge \text{Cov}(z, y)\text{Var}(y)^{-1} \text{Cov}(y, z)=\left[\frac{\partial}{\partial \theta} \psi(\theta)\right]^TI(\theta)^{-1}\frac{\partial}{\partial \theta} \psi(\theta)$$
