I have written some code that can do Kalman filtering (using a number of different Kalman-type filters [Information Filter et al.]) for Linear Gaussian State Space Analysis for an n-dimensional state vector.
Let us consider the linear Gaussian State Space model
$$y_t = \mathbf{Z}_{t}\alpha_{t} + \epsilon_{t},$$ $$\alpha_{t + 1} = \mathbf{T}_{t}\alpha_{t} + \mathbf{R}_{t}\eta_{t},$$
where $y_{t}$ is our observation vector, $\alpha_{t}$ our state vector at time step $t$. The quantities in bold are the transformation matrices of the state space model which are set according to the characteristics of the system under consideration. We also have
$$\epsilon_{t} \sim NID(0, \mathbf{H}_{t}),$$ $$\eta_{t} \sim NID(0, \mathbf{Q}_{t}),$$ $$\alpha_{1} \sim NID(a_{1}, \mathbf{P}_{1}).$$
where $t = 1,\ldots, n$. Now, I have derived and implemented the recursion for the Kalman Filter for this generic state space model by guessing the initial parameters and variance matrices $\mathbf{H}_{1}$ and $\mathbf{Q}_{1}$.
where the points are the Nile River water levels for Jan over 100 years, the line is the Kalamn Estimated state, and the dashed lines are the 90% confidence levels.
Now, in order to estimate the initial parameters of my model ($\mathbf{H}_{1}$ and $\mathbf{Q}_{1}$) I have adopted the Expectation Maximisation method to maximise the loglikelihood function
$$\log L(Y_{n}|\psi) = -\frac{np}{2}\log(2\pi) - \frac{1}{2}\sum^{n}_{t = 1}(log|\mathbf{F}_{t}| + v^{\mathsf{T}}_{t}\mathbf{F}_{t}^{-1}v_{t}), \;\;\;\; (1)$$
with respect to parameter vector $\psi$. Where $v_{t}$ is the state error and $\mathbf{F}_{t}$ is the state error variance. The vector $\psi$ contains the unknown parameters (i.e., either the measurement and state disturbance variances themselves, or a re-parametrisation of these same variances, see below)
A crucial role in the maximisation of (1) is played by the gradient or score vector
$$\partial_{1}(\psi) = \partial \log L(Y_{n}|\psi)/\partial\psi,\;\;\;\; (2)$$
At least two algorithms can be used to maximise (1). The first is the so-called EM (Expectation-Maximisation) algorithm, and the second is the BFGS (Broyden-Fletcher-Goldfarb-Shannon) algorithm. I have successfully implemented the EM algorithm and this is providing correct results but is very slow; in light of this I now wish to implement the BFGS algorithm as this has better performance.
So, from (2) using Taylor's theorem we can write
$$\partial_{1}(\psi) \simeq \tilde{\partial}_{1}(\psi) + \tilde{\partial}_{2}(\psi)(\psi - \tilde{\psi}),$$
where $\partial_{2}(\psi)$ is the Hessian matrix or second derivative of the loglikelihood function wrt $\psi$ and $\tilde{\partial}_{N}(\psi) = \partial_{N}(\psi)|_{\psi = \tilde{\psi}}$. So now we can write our updated parameter vector as
$$\bar{\psi} = \tilde{\psi} - \tilde{\partial}_{2}(\psi)^{-1}\tilde{\partial}_{1}(\psi) = \tilde{\psi} - s\tilde{\pi(\psi)}.\;\;\;\; (3)$$
where here $\tilde{\pi(\psi)} = - \tilde{\partial}_{2}(\psi)^{-1}\tilde{\partial}_{1}(\psi)$ and $s$ is our step size. So the basic outline of the BFGS algorithm to get our parameter vector and maximize (1) using BFGS is:
Initialise parameter vector $\psi = \psi^{*}$.
Then apply the Kalman and disturbance smoothing filters and thus for score vector (2) at $\psi = \psi^{*}$. When running the Kalman Filter evaluate the log-likelihood function.
Use (3) to obtain new values for $\psi$ given by $\psi^{+}$ and go back to step 2 until the value of loglikelihood function (1) no longer improves.
This seems fine, but my problem is how to calculate the step size $s$ in (3) above. I know a line search is required, but remember $\psi$ is a vector and it is not clear what function I can use in a traditional line search technique in order to evaluate $s$.
Question: what line search method should I use in order to evaluate $s$?
Ps. Assume that I would like to use the zoom line search algorithm (algorithm 3.5 of Numerical Optimization by Nocedal & Wright).
Further to the very helpful answer by @LKlevin, I want to classify my problems more explicitly...
$^{*}$Most (L-)BFGS methods accept vector functions, that is a function $f$ that depends on a n-dimensional (in this case let us assume two-dimensions) parameter vector $\mathbf{x} = (x_{1}, x_{2})^{\mathsf{T}}$, so if we had $F(x_{1}, x_{2}) = 2x_{2}^{2} + 4x_{2}$ you would pass into the BFGS algo. $F$ and $\nabla{F}$. This is fine and the text I reference as a PDF above tells me all I need to know here in terms of using BFGS to minimise this function.
However, in the case I have above, from (3) we have
$$\bar{\psi} = \tilde{\psi} - \tilde{\partial}_{2}(\psi)^{-1}\tilde{\partial}_{1}(\psi) = \tilde{\psi} - s\tilde{\pi(\psi)},$$
which can be written as
$$\left[\begin{array}{c} \psi_{1}\\ \psi_{2}\\ \end{array}\right] = \left[\begin{array}{c} \tilde{\psi}_{1}\\ \tilde{\psi}_{2}\\ \end{array}\right] - \mathcal{H}\,\partial_{1}(\psi)\left[\begin{array}{c} \tilde{\psi}_{1}\\ \tilde{\psi}_{2}\\ \end{array}\right] = \left[\begin{array}{c} \tilde{\psi}_{1}\\ \tilde{\psi}_{2}\\ \end{array}\right] - s\left[\begin{array}{c} \tilde{\pi}_{1}(\psi_{1})\\ \tilde{\pi}_{2}(\psi_{2})\\ \end{array}\right],$$
where $\mathcal{H}$ is the inverse Hessian matrix.
So the problem seems to consist of $n$ linear optimizations, not a single equation depending on a vector of parameters. My problem here is the use of the vectors $\psi$ in this case above.
So, can I optimise this without splitting the n-dimensional (2 dimensional in the example above) vector into $n$ equations which I can optimise using functions as described in the paragraph $^{*}$ above?
Another issue is that the loglikelihood function is not explicitly dependent on $\psi$, so how can i use this in the line search method?
