Does Fisher scoring exist as such? I'm studying second order optimization methods in statistics, and I've run into a conceptual barrier that I was hoping someone can help me with.
First for some notation: consider a sample $X_1,\dots,X_n\thicksim P_{\theta^\star}$ [edit: $\theta^\star$ is the true, unknown parameter. I had originally denoted it as $\theta_0$, which has been misunderstood as the first of the iterates $\theta_k$ discussed below -- sorry for the poor notation], and our goal is to find the MLE $$\hat\theta=\operatorname{arg\,max}_{\theta}L(\theta)=\operatorname{arg\,max}_{\theta}\sum_{i=1}^n\ell(X_i, \theta)=\operatorname{arg\,max}_{\theta}\sum_{i=1}^n\log P_{\theta}(X_i)$$.
We consider two iterative algorithms:

*

*Newton-Raphson: $\theta_{k+1}=\theta_k + (-\nabla^2 L(\theta_k))^{-1}\nabla L(\theta_k)$

*Fisher scoring/natural gradient: $\theta_{k+1}=\theta_k + (I(\theta_k))^{-1}\nabla L(\theta_k)$, where under the usual assumptions $$I(\theta)=E_{\theta^\star}[-\nabla^2 L(\theta)].$$
[edit:the expectation in this formula must be over $\theta$ instead of over $\theta^\star$, this misunderstanding was at the core of my confusion.]
Now, it seems to me that: since we don't know $\theta^\star$ (and therefore expectations under $P_{\theta^\star}$ can't be computed), it follows that to implement Fisher scoring we must approximate $I(\theta_k)$ with some approximation, for instance $(\nabla L(\theta_k))^\top(\nabla L(\theta_k))$ or $-\nabla^2 L(\theta_k)$. The first approximation is reminiscent of the Gauss-Newton algorithm, and the second is exactly Newton-Raphson.
Assuming what I have written above is correct/makes sense, my question is: am I correct that in implementation, Fisher scoring is always carried out by an approximation like the above? Is it ever possible to carry out "exact" Fisher scoring using the expected information matrix? I've read a lot of sources where these algorithms are presented separately, so I'm uncertain, and I'd appreciate any clarification. Sorry if this is a duplicate -- there are lots of similar questions here which I have read, but I could not find one that addressed my question.

Yes, certainly it exists and there is no problem with its implementation.
Fisher scoring is one of the most commonly-used algorithms in statistics and is it usually implemented exactly as defined in your equation. It is by far the most frequently used algorithm for generalized linear models or for non-linear least squares, for example.
Iterative algorithms like Fisher scoring and Newton-Raphson need to be started from an initial parameter estimate. Any starting value can be used, and one of the advantages of Fisher scoring is that it tends to converge to the correct maximum likelihood estimate from a wide range of starting values. The starting value has to be chosen separately from the iterative algorithm. Typically a moment estimator might be used. For generalized linear models, Fisher scoring is usually started by writing the interation in terms of fitted values and substituting the observed responses in place of the initial fitted values---that then generates $\theta_0$ from linear regression of the link-transformed responses on the covariates. In the worst case, you could simply try starting the algorithm from an arbitrary value or from a randomly chosen value. Whatever the starting value is, it is represented by $\theta_0$ in the definition of the iteration. In the algorithm definition, $\theta_0$ simply represents any given value that belongs to the interior of the parameter space.
In Fisher scoring, the information matrix is computed analytically and exactly. The information matrix, either observed (for Newton-Raphson) or expected (for Fisher scoring) is generally computed exactly at the current working value of $\theta_k$.
In the context of nonlinear regression with normally distributed responses, Fisher scoring and Gauss-Newton are identical.
Further remarks

*

*Some remarks on notation. The true parameter value is simply written as $\theta$. The iteration is started at $k=0$ so $\theta_0$ is the starting parameter value. Assuming the iteration converges, the maximum likelihood estimator is $\hat\theta=\theta_\infty$. Please don't try to use $\theta_0$ for the true value, that is a serious abuse of notation in this context that will just confuse everyone!


*We do not have an iid sample $X_1,\ldots,X_n$. Fisher scoring is always applied in regression contexts in which every $X_i$ has a different distribution. In other words, we have a single vector random variable $X$ of dimension $n$, and $\theta$ is a vector parameter than applies to $X$ as a whole. You cannot write down the log-likelihood function as a sum of terms for individual $X_i$ unless you make reference to the different covariate values that apply to each $X_i$.


*Likelihood derivatives are computed with respect to the log-likelihood $\ell$, not with respect to the unlogged likelihood $L$.


*The formula for $I(\theta)$ is established assuming that $\theta$ is the true value, i.e., we evaluate $I(\theta)= E_\theta(-\nabla^2\ell(\theta)$ assuming that $\theta$ is the true value. (Note that the three values of $\theta$ in this formula are always the same, one cannot mix three different $\theta$s as you do in your question!) This establishes $I()$ as an analytic mathemical function that can then be evaluated at any other value of $\theta$. When we evaluate $I(\theta_k$) in the Fisher scoring algorithm, we are not reevaluating the likelihood expectation or assuming that $\theta_k$ is the true value, we are just applying $I()$ as a mathematical function. When applying the algorithm, there is no assumption that $I(\theta_k)$ is the true expected information, it is just a numerical formula. It is all much simpler than you seem to be assuming.
References
Smyth, G. K. (1998). Optimization and nonlinear equations. In: Encyclopedia of Biostatistics, P. Armitage and T. Colton (eds.), Wiley, London, pp. 3174-3180. http://www.statsci.org/smyth/pubs/EoB/bao021-.pdf
Smyth, G. K. (2015). Optimization and nonlinear equations. Wiley StatsRef: Statistics Reference Online (edited by W. Piegorsch), John Wiley & Sons, https://doi.org/10.1002/9781118445112.stat05030.pub2, pages 1-9. http://www.statsci.org/smyth/pubs/OptimNonlinEqnPreprint.pdf
