My question is how to solve specified matrix equation (see bellow). However let me first explain background and where the equation comes from.

Kalman filter allows us to estimate state at time $t$ as multivariate normal distribution and the formulas for kalman filter (see wikipedia) says what is the mean and the covariance of this distribution.

The covariance matrix of estimated state doesn't depend on values of observations (but depends on their covariance matrix). The situation is even simper when this covariance and any other the parameters of Kalman filter doesn't change over time.

My question is: What is limit covariance of state estimate after infinite time has passed?

This question is equivalent of asking what is the limit of $\lim_{x\to\infty} P_t$ where

$$P'_t = F P_{t-1} F^T + Q$$

$$P_t = (I-P'_t H^T (R+H P'_t H^T)^{-1} H) P'_t$$

$I$ is unit matrix.

Under some conditions the question is equivalent to solving the above equation with $P_t = P_{t-1}$.

I believe that the series converges under reasonable assumtions to the solution of the equation.

Note that all $P_t$, $Q$ and $R$ are symetric possitive definite and all eigenvalues of F are less or equal to one.

How do I solve the equation or find limit of $P_t$?

• Stefan Kohl: I am new stacks exchange. What do you mean by top level tag? Where should I add it? What should I do? Sep 18 '17 at 17:39
• Aha, maybe I understand. I changed the tags. Sep 18 '17 at 18:31
• Nothing commutes in general, right? If so, I'm extremely skeptical about the possibility to write a neat closed form expression for the solution but we can still try to investigate the convergence and figure out some other things you might be interested in. So, can you tell us what you'll be happy with in lieu of an explicit formula? Nov 25 '17 at 3:35

Under the assumption that the pair $$(F,H)$$ is detectable then there is a positive definite solution where $$P=P_t=P_{t-1}$$. I assume that with "all eigenvalues of F are less or equal to one" you are referring to the absolute value of the eigenvalues. In that case $$F$$ would be a Schur matrix and thus $$(F,H)$$ would always be detectable for any $$H$$.

Solving $$P_t=P_{t-1}$$ can be done using any discrete algebraic Riccati equation (DARE) solver. Your stationary equation can be formulated into the following more standard representation of a DARE

$$A^\top P\,A - P - A^\top P\,B \left(B^\top P\,B + \bar{R}\right)^{-1} B^\top P\,A + \bar{Q} = 0,$$

with

$$\bar{R} = R + H\,Q\,H^\top,$$

$$\bar{Q} = Q - Q\,H^\top \bar{R}^{-1} H\,Q,$$

$$A = F^\top - F^\top H^\top \bar{R}^{-1} H\,Q,$$

$$B = F^\top H^\top.$$

A common algorithm used to solve a DARE is discussed in the paper: Laub, Alan. "A Schur method for solving algebraic Riccati equations." IEEE Transactions on automatic control 24, no. 6 (1979): 913-921 PDF.

One way to make this problem more tractable is to examine the limit of continuous time instead of discrete time. The scaling of the transformations in the continuum limit should be

$$F = 1+fdt\\ Q = qdt\\ R = r/dt\\ H = h$$

Plugging those in to the steady state equations and keeping terms up to first order in $$dt$$, you should get

$$fP+Pf^T+q-Ph^Tr^{-1}hP=0,$$

which is hopefully easier to work with.

The special case $$f=0, h=1$$ gives the matrix geometric mean

$$P = q^{1/2}(q^{-1/2} r q^{-1/2})^{1/2}q^{1/2}\text,$$

and it looks like you can get the general case into that form too by completing the square.