Suppose $l(b; y)$ is the log-likelihood function we want to maximize, i.e. we want to solve $U = l' = 0$ for the ML estimation of $b$. Using Newton-Raphson method, we have
\begin{align*} b^{(m)} = b^{(m-1)} - \frac{U}{U'}, \end{align*}
This can be again approximated by: \begin{align*} b^{(m)} = b^{(m-1)} - \frac{U}{E(U')}, \end{align*}
If the random variable $y$ is from a distribution belonging to the exponential family, then we have $M = Var(U) = -E(U')$, thus:
\begin{align*} b^{(m)} = b^{(m-1)} + \frac{U}{M}. \end{align*}
If $\mathbf{b}$ is a vector of several parameters, and $\mathbf{U}$ is vector of derivatives where $U_j = \frac{dl}{db_j}$, then it's said the above formula can be generalized to
\begin{align*} \mathbf{b}^{(m)} = \mathbf{b}^{(m-1)} + \mathbf{M}^{-1} \mathbf{U}, \end{align*}
Where $\mathbf{M}$ is the variance-covariance matrix of $\mathbf{U}$
I am wondering how this generalization can be derived, book recommendations are also very appreciated.