Extended Kalman filtering on Stiefel manifolds

1 Introduction

There are multiple situations where projected dynamics is inherently tied to gaining information. We shall begin by giving simple examples of two such situations in order to broaden the perspective before we go into the specifics relating to applications of Stiefel manifold-valued statistics in radiology.

Consider a radioactive dust particle travelling through an open space. Suppose also there is a moderate wind going through the area and we can only measure the direction of the incoming radiation. In physics one often models the movement of dust particles as if they travel like Brownian motion with a drift (the drift here comes from the wind). Therefore the dust particle does not travel deterministically and there is an uncertainty inherent to the position of the dust particle. Moreover, the measured radiation incoming into our sensor might have a somewhat perturbed inclination error as it collides with air particles in between the dust particle and the sensor. Hence, there is an uncertainty of the incoming direction of the radiated particles as well. In this situation the dust particle is travelling across R 3 and our measurements are taken on the sphere S 2. The information obtained is (non-linearly) projected onto the sphere.

Another plausible scenario would be a circular particle accelerator. As a charged particle is centripetally accelerated through electro-magnetic forces we measure its position in the accelerator using sensors throughout the walls of the accelerator. As the particle travelling through the accelerator is a quantum particle, the position of the particle is inherently probabilistic. At the same time the sensors are not perfect and there is uncertainty in the measurements. Here the particle moves inside a tubular neighbourhood of a circle and the measurements happen on the torus T 2.

In MIMO radio-systems for cellular networks a configuration of n antennas interfere constructively and a directed radio signal is produced with high energy efficiency. This signal is then received at and responded by k receivers. These receivers move around, in part, in a predictive way (a person moving in a certain direction will probably move in that direction for a while) and, in part, with some uncertainty. Hence, these receivers may be modelled as stochastic processes with some dynamical component. Moreover, the antenna configurations can be modelled as being a vector in R n for each receiver. In total we consider the state as an n by k matrix. However, as these configurations do not depend on the total strength of the currents, only their relative strength and phase, so measurements are realised as n by k orthogonal matrices, see Hussien et al. (2014); Pitaval and Tirkkonen (2013); Seddik et al. (2017). These realisations may be computed by only retaining the orthogonal part, Q, in the polar decomposition of the n by k matrix. This Q can be computed so that it is a projection onto the space of n by k orthogonal matrices, i.e. the Stiefel manifold S t n , k. These measurements are assumed to be noisy and the measured optimal configuration of the antennae has some uncertainty.

All these three situations illuminates a mathematical problem: an object moves around with a drift in some (possibly an open subset of) vector space, see Figure 1 for a picture of these situations. If one performs measurements on a non-linear space, i.e. if x , y ∈ X, then it is generally not true that x + y ∈ X. How would one filter out this information using both the known drift of the system and the measurements? In all the above examples the non-linear spaces are manifolds which have sufficient geometrical structure which allows for measurement of distances between points. Using distances and curves on these manifolds it is then possible to weight the prediction together with the measurement. In a nutshell, working on the original filtering equations developed in Kálmán (1960a); Kálmán (1960b); Kálmán and Bucy (1961) we will filter the predicted point and the measurement by using a curve that connects these two points. Then one can obtain the filtered mean by taking a weighted average along the curve.

In recent time, quite a few manifold-valued filters has been developed. Some Snoussi and Mohammad-Djafari (2006); Chahbazian et al. (2022); Boulanger et al. (2016); de Figueredo et al. (2022) apply some type of particle filters in order to solve the manifold-valued filtering problem. However, particle-filters are accurate but computationally expensive, especially as the dimension of the manifold goes up. For example, problems in radiology where real time estimations of (ideally ≥ 100 dimensional) manifolds the particle filters are non-tractable. In Boulanger et al. (2016) the authors explicitly explain the curse of dimensionality of these methods by saying

“For higher dimensions ( n > 10 ), the use of a particle filter should not be considered. Indeed, as n increases, the dimension of the estimation space increases like O ( n 2 ) and a very large number of particles is required to perform a proper estimation”.

Any Manifold-valued Extended Kalman filters that has been previously considered Brossard et al. (2021); Barrau and Bonnabel (2018); Barrau and Bonnabel (2016) only consider the case when the manifold itself is a Lie group or by lifting the problem to a Lie group Mahony et al. (2022). Therefore, there is a need to develop some kind of Extended Kalman filter with projected normal observations on the Stiefel manifolds that can also deal with high dimensions.

More broadly, this paper can be viewed as part of the vast field of directional statistics. The area of directional statistics has found applications in medicine, see Craig and Oppe (2010); Demir and Bilgin (2019); Karaibrahimoglu et al. (2021), in meteorology, see Li et al. (2018); Nuñez-Antonio et al. (2015); Nuñez-Antonio et al. (2011), and in robotics, see Lang and Feiten (2012); Sveier and Egeland (2018), to name a few. Applications of these manifold-valued filtering methods are found in diverse fields. In aerospace and navigation, they are utilized for spacecraft attitude stabilization and control Bloch and Marsden (1990); Astolfi and Lovera (2002) as well as terrain-based positioning Murangira et al. (2016). Robotics and sensor networks employ these techniques for rigid body localization Chepuri et al. (2014) and kinematic state estimation Park et al. (2008). In computer vision and medical imaging, they facilitate action and object recognition Turaga et al. (2011); Lui (2012); Liu et al. (2004) and the statistical analysis of diffusion tensor MRI Lenglet et al. (2006). Furthermore, specific applications are found in biomechanics for assessing climbing fluency Seifert et al. (2014), in wireless communications for predictive quantization Schwarz and Rupp (2014), and in deep learning for enforcing orthogonality constraints in neural network training Bansal et al. (2018); Wang et al. (2020); Fiori (2001); Fiori (2008).

This paper is mainly aimed towards the engineering/application viewpoint of the Extended Kalman filter with Stiefel manifold-valued observations. While detailed convergence proofs are provided in a sister paper Figueras et al. (2025a), we present the mathematical derivation and Bayesian motivation for the filter design in Section 3. However, in order to make any technical conclusions one has to make a few approximations of some of the maps used in the filter. These approximations are not necessary in order to apply our proposed filter and the approximations give rise to additional error propagation. Therefore, we keep the more exact formulation and application relevant version of the filter in this paper with emphasis on application and proof of concept.

The paper is organized as follows. Section 2 explains the foundational concepts for the problem at hand. In Section 3 we give an exposition on how to use the extended Kalman filter when the measurements are points on the Stiefel manifold S t n , k. The reader with little interest of the geometrical meaning behind the Extended Kalman filter on S t n , k may jump directly to Algorithm 2 at the end of this section. In Sections 4 and 5 we give results from simulations for the 2-sphere and S t 4,2, respectively. In Section 7 we discuss some limitations and we present some potential drawbacks of the extended Kalman filter on Stiefel manifolds.

2 Setting

Suppose an object X t is represented as an n by k matrix, that is X t ∈ R n × k. As t increases we suppose it moves around with known dynamics and with non-zero noise (therefore the movement is not fully deterministic). This may be mathematically represented as an SDE (Stochastic Differential Equation)

d X t = A X t d t + ν d B t , X 0 = d N μ 0 , ν 0 2 i d R n × k . (1)

Here i d R n × k denotes the identity map from the space of n by k matrices to itself, i.e. for given Y ∈ R n × k, i d R n × k ( Y ) = Y. More explicitly, if Y vec ∈ R n k is the vectorization of an n times k matrix Y, then if i d R n × k vec is the corresponding vectorized linear map, it is the n k × n k identity matrix. The initial value X 0 is the starting (normal) distribution of X t and A ∈ R n × n describes the dynamics which is the infinitesimal displacement over time. The term d B t is a formal stochastic differential of the Brownian motion which we scale by a factor ν. After a time step δ t > 0 and if A = 0, then solving Equation 1 one would obtain X t = X 0 + ν ∫ 0 δ t d B t which has variance matrix ( ν 0 2 + ν 2 δ t ) i d R n × k and constant mean μ 0. More generally, if A is anti-symmetric, i.e. A T = − A, then one obtains the full solution to Equation 1 asHere exp M : R n × n → R n × n denotes the matrix exponential defined aswhere A j denotes j-times matrix multiplication and A 0 = I n × n, the n × n identity matrix. (Recall that this series is absolutely convergent for any A).

X t = d N exp M t A μ 0 , ν 0 2 + t ν 2 id R n × k .

exp M A = ∑ j = 0 ∞ A j j ! ,

Remark 2.1 Note that this solution is only true when the factor ν is a scalar (for which we then call the noise isotropic). If the factor ν is not a scalar then there will appear more time dependent terms in the covariance part.

Suppose now the process X t from Equation 1 is observed, but we only observe the orthogonal part in the polar decomposition of X t. This observation is a projection onto the space of n × k matrices with pairwise orthonormal columns, the Stiefel manifold

S t n , k ≔ X ∈ R n × k : X T X = I k .

We write p r : R n × k → St n , k, and its computation can be done by a polar decomposition using SVD algorithms. If X t = Q t D t S t * is the thin SVD-decomposition, see Bai et al. (2000), where D t and S t are k × k matrices and Q t is n × k, then p r ( X t ) = Q t S t *. Note that this map is heavily non-linear and should not be confused with (linear) projections seen in linear algebra. A conceptual picture on how this projection of the stochastic process X t might look like is shown in Figure 2.

Suppose further that when this projected process is observed there is an error of magnitude ξ 2 inherent to these measurements. That is, if p r ( X t ) is the orthogonal part of X t using the polar decomposition, the measurement’s error of X t is assumed to have mean p r ( X t ) and (scalar) variance ξ 2 (we shall delve into what this means in the next section). In Figure 1 there is an illustration of how these faulty measurements may look like.

Now the goal is to estimate p r ( X t ) given that we know the system parameters A ∈ R n × n, ν ∈ R and the measurement error. In order to capture parts of this non-linearity we shall use the tangent space. The tangent space T X S t n , k contains all possible velocities (aka. tangent vectors) at the point X ∈ S t n , k such that any curve going through X with velocity in T X S t n , k will (infinitesimally) stay on S t n , k, see Figure 3 for an illustration. Explicitly, the tangent space is given bywhich may be deduced from that any curve Y ( t ) on S t n , k must satisfy Y T ( t ) Y ( t ) = I k. Differentiating this expression yields exactly the above set of tangent vectors.

T X S t n , k = V ∈ R n × k : V T X = − X T V

Now, considering S t n , k as a subset ¹ of R n × k, one can measure distances on S t n , k by considering the length of curves on S t n , k with constant velocity by restricting a modified version of the Frobenius norm onto the tangent space at every point. Recall that the classical Frobenius norm of a matrix B ∈ R n × k is defined by ∥ B ∥ 2 ≔ ‖ B ‖ F 2 = t r ( B T B ). Through the notion of distance, one can relate points on the tangent space T X S t n , k to points near X: by denoting exp X ( V ) as the point Y ∈ S t n , k which is reached by taking the constant velocity curve starting at X with starting velocity direction V / ∥ V ∥ and is the distance ∥ V ∥ away. Note here that exp X is not the same as the matrix exponential map exp M defined earlier. The map exp X : T X St n , k → St n , k is called the (Riemannian) exponential map (we have provided a picture of this exponential map in Figure 4). We denote the inverse of exp X by log X, and call it the logarithm map.

Remark 2.2 Explicitly, one may compute exp X ( V ) by first computing the QR-decomposition (here Q is a n by k orthogonal matrix and R is a k by k upper triangular matrix), then

In Theorem 2.1 and Corollary 2.2 in Edelman et al. (1998) one can find more information about this exponential map. To compute the logarithm, which is a little bit more involved to compute but still possible, we refer the reader to Zimmermann (2017).

The inner product, aka metric, on the tangent space T X S t n , k, with X ∈ S t n , k, is computed as follows

⟨ V , W ⟩ X = t r V T I n × n − 1 2 X X T W ) . (2)

Remark 2.3 An important key point to note is that S t n , k is not considered as a subset of R n × k with the standard Frobenius norm. Here, the norm induced by Equation 2 is the norm inherited from the Frobenius norm on the Lie group O ( n ) of which the Stiefel manifold is a quotient space. This actually makes S t n , k a isometrically submerged submanifold of O ( n ) and by extension R n × n. It is still an open question how and if one can isometrically embed S t n , k in R n × k with an easily expressed norm defined on R n × k.

Through the logarithm map one can define a notion of average which is inherent (aka. intrinsic) to S t n , k. Given points Y 1 , … , Y N ∈ S t n , k, their average is defined as Y ̄ ∈ S t n , k which minimises the distance squared from all the points Y 1 , … , Y N. Locally, the average Y ̄ ∈ S t n , k explicitly satisfies the barycenter equation

∑ i = 1 N log Y ̄ Y i = 0 .

This average is illustrated in Figure 5. Note that this is a non-linear generalisation of the standard definition of average in a vector space: in this case the average y ̄ ∈ R satisfieswhich is a bit more convoluted way of saying

1 N ∑ i = 1 N y i − y ̄ = 0

y ̄ = 1 N ∑ i = 1 N y i .

Let B j be an orthonormal basis for the vector space T Y ̄ S t n , k under the inner product Equation 2. The sample covariance given points Y 1 , … , Y N ∈ S t n , k, is then the n k − 1 2 k ( k + 1 ) by n k − 1 2 k ( k + 1 ) matrixwith corresponding scalar variance

C o v Y i , j = 1 N − 1 ∑ ℓ = 1 N 〈 log Y ̄ Y ℓ , B i 〉 Y ̄ 〈 log Y ̄ Y ℓ , B j 〉 Y ̄

t r C o v Y n k − 1 2 k k + 1 . (3)

Note also that the notion of variance has a meaning of measuring the squared distance, and the mean has the meaning of being the point which minimises the least squared distance, see (Figueras et al., 2025b, Definition 2.1 and 2.2).

Observe that spheres S n − 1 are special cases of S t n , k. Specifically note that S n − 1 is the same as S t n , 1. On the spheres, we know that given a normal random variable X in R n + 1 with mean μ and isotropic covariance σ 2 I n + 1 the projected random variable p r ( X ) = X ∥ X ∥ has inherent average μ ‖ μ ‖. Moreover, it has been shown that there is unique, one-to-one, intrinsic scalar variance of p r ( X ) for every choice of σ for S 2, see (Figueras et al., 2025b, Theorem 2.8). Along this paper we shall assume that the same is true for all S t n , k.

Let X ∈ R n × k is a normal random variable with mean E [ X ] ∈ S t n , k and covariance matrix C o v ( X ) = σ 2 id R n × k. Then the inherent scalar variance can be computed bywhich is the true scalar variance of the random variable p r ( X ) corresponding to the sample scalar variance in Equation 3. Because the intrinsic covariance is inherently respecting the curvature of the manifold, the mapping η gives a (virtually) cost-free way of taking the curvature of the manifold into consideration when moment-matching the covariance. This allows to translate between true scalar variance and the corresponding variance in the ambient space, see Section 3. Due to Figueras et al. (2025b), in the case of S t 3,1 = S 2, η is known to be injective and is explicitly given by

η σ 2 = 1 n k − 1 2 k k + 1 ∫ St n , k dist 2 E X , x p pr X x dVol St n , k x , (4)

η σ 2 ≔ 1 2 π 1 / 2 exp − 1 2 σ 2 ∫ 0 π ϕ 2 cos ϕ σ + cos 2 ϕ σ 2 + 1 Φ cos ϕ σ φ cos ϕ σ sin ϕ d ϕ . (5)

As of today, there is no such explicit formula for a general Stiefel manifold since p p r ( X ) is difficult to express for S t n , k when k ≥ 2. Because of this, we shall need the following conjecture in order to perform parameter inference on the projected observations.

Conjecture 2.4 Let X ∈ R n × k be a normal random variable with E [ X ] ∈ S t n , k and covariance C o v ( X ) = σ 2 id R n × k . Then the projected normal random variable has inherent average E [ X ] . Moreover, there exists an injective function η : [ 0 , ∞ ) → [ 0 , ∞ ) such that the inherent scalar variance is uniquely determined by η ( σ 2 ) for each σ > 0.

Remark 2.5 Moreover, the function η is conjectured to have derivative 1 at 0. This has been numerically verified as seen in Figure 6. This is also a reasonable conjecture intuitively. When the variance is small then the projected normal lives in a small neighborhood on the manifold and the manifold looks flat.

The pseudo-code in Algorithm 1 help us to numerically estimate the function η using the Monte Carlo method, in the case that η is difficult to explicitly express. That Algorithm 1 converges follows from (Figueras et al., 2025b, Theorem 2.5).

Algorithm 1 An outline of how to estimate the function η, defined by Equation 4, on St_*n,k*.

Using Algorithm 1, we have numerically verified Conjecture 2.4 for the Stiefel manifolds S 2, S t 6,2, S t 8,2 and S t 8,3, see Figure 6. For each σ 2 ∈ ( 0,1 ] and each Stiefel manifold, 100 pseudo-random samples of N ( I n , k , σ 2 i d R n × k ) were drawn. Then, η ( σ 2 ) was estimated for each σ 2 and each Stiefel manifold. Notice that S 2 is the only exception, in this case η ( σ 2 ) was computed by numerically approximating the integral in Equation 5. Notice also that each η is roughly monotonically increasing and does not contradict Conjecture 2.4.

Now we have sufficient amount of tools to do filtering on S t n , k.

3 Extended Kalman filtering on Stiefel manifolds

Suppose we have the following filtering problemwhere t m is an increasing sequences of discrete times, ε m ∈ T p r ( X t m ) S t n , k are i. i.d. such that ε m = d N ( 0 , ξ 2 i d T pr ( X t m ) St n , k ), ν ∈ R and A : R n → R n is an anti-symmetric linear mapping ( A T = − A ).

d X t = A X t d t + ν d B t , X 0 = d N μ 0 , σ 0 2 i d R n × k Z m = p r pr X t m + ε m , (6)

Next we shall adopt the classical Extended Kalman filter to the setting of Stiefel manifolds. Suppose we have prior X 0 = d N ( μ 0 , σ 0 2 i d R n × k ), and suppose after time t 1 we make measurement z 1. Now our predicted average is μ pred. = exp M ( t 1 A ) μ 0 and the (inherent) predicted scalar variance is, by a linear approximation, P pred. = η ( σ 0 2 + t 1 ν 2 ). Since the measurements have scalar variance ξ 2 we can compute the Kalman gain by

K = σ 0 2 + t 1 ν 2 σ 0 2 + t 1 ν 2 + ξ 2 .

We use the logarithm map to get the tangent vector y 1 ∈ T μ pred. S t n , k which gives the velocity of the constant speed curve starting at μ pred. which reaches z 1, i.e. exp μ pred. ( y 1 ) = z 1. Now we rescale this vector by the Kalman gain K to get a weighted average of the predicted point μ pred. and the measurement z 1. We end up with the mean estimate exp μ pred. ( K y 1 ) and the estimate ( 1 − K ) P pred. for the inherent (scalar) variance. Hence our posterior, after lifting this distribution back to the ambient space R n × k, is

N exp μ pred. K y 1 , η − 1 1 − K P pred. id R n × k .

In Figure 7 we give a picture of the geometry behind the filtering. Moreover, the entire extended Kalman filtering algorithm is given step-by-step in Algorithm 2 as pseudo-code.

3.1 Mathematical interpretation

The proposed filter can be interpreted within the framework of Bayesian estimation on manifolds. The filtering problem aims to recursively compute the posterior probability distribution of the state X t given the sequence of measurements z 1 , … , z m. Due to the non-linear nature of the projection map p r ( ⋅ ) and the Stiefel manifold geometry, the exact posterior is generally non-Gaussian and computationally intractable to represent exactly.

Our approach approximates the posterior distribution at each step k by a Projected Normal distribution PrN ( μ k , P k ), which is characterized by a central point μ k ∈ S t n , k and a covariance matrix P k defined on the tangent space T μ k S t n , k.

3.2 Bayesian interpretation

The filter operates in a predict-update cycle analogous to the standard EKF:

• Prediction: The predictive distribution is obtained by propagating the previous estimate through the system dynamics. The mean is updated via the noiseless dynamics μ pred = exp M ( t A ) μ k − 1, and the covariance is inflated to account for process noise, approximated here by the function η which maps ambient isotropic variance to intrinsic manifold variance.
• Update: The update step assimilates the new measurement z m. This is performed in the tangent space of the predicted mean μ pred. The measurement z m is mapped to the tangent vector y = log μ pred ( z m ). The filter then computes the maximum a posteriori (MAP) estimate of the linearized model in the tangent space. This results in a Gaussian posterior on the tangent space N ( K y , ( 1 − K ) P pred ), which is then retracted back to the manifold to form the new state estimate μ K = exp μ pred ( K y ).

3.3 Error criterion

The update step explicitly minimizes the mean squared error (MSE) of the linearized error in the tangent space. Specifically, the Kalman gain K is derived to minimize the trace of the posterior covariance matrix in the tangent space. This corresponds to minimizing the quadratic cost function:where v ∈ T μ pred S t n , k represents the correction to the predicted mean, and R = ξ 2 i d is the measurement noise covariance. By minimizing the cost defined by Equation 7, the filter finds the optimal trade-off between the prior information (prediction) and the new information (measurement) consistent with the local Riemannian metric.

J v = v T P pred − 1 v + y − v T R − 1 y − v , (7)

Remark 3.1 The update step in Algorithm 2 is theoretically justified by interpreting the filter as a recursive maximum a posteriori (MAP) estimator acting on the tangent space T μ pred. S t n , k. By defining a quadratic cost function J ( v ) that penalizes weighted deviations from both the predicted prior and the measurement innovation, the optimal state correction v * is found at the point where the gradient of the log-posterior vanishes. Applying the Matrix Inversion Lemma to the resulting linear system transforms the solution into the standard Kalman gain form v * = K y, where K = σ 0 2 + t ν 2 σ 0 2 + t ν 2 + ξ 2, confirming that the retraction of the weighted tangent vector back to the manifold provides the optimal state estimate within the local Gaussian approximation.

Algorithm 2 One step of the extended Kalman filter on a manifold after time t.

4 Simulations on the 2-sphere

Throughout we consider system parameters ν 2 = 1, and ξ 2 = 0.1 with initial covariance Σ 0 = 0.1 I 3 and initial mean

A = 0 0.263 0.036 − 0.263 0 − 0.653 − 0.036 0.653 0 , (8)

μ 0 = 0 0 1 . (9)

Now consider stochastic process X t solving Equation 1. The system process X t is simulated as follows: Consider constant interval discrete times 0 = t 0 , … t N where Δ t = t j − t j − 1 for 1 ≤ j ≤ N. We choose N = 2000. The initial point is drawn from X 0 = N ( μ 0 , Σ 0 ). Now each point X t j is simulated by drawing V j from N ( 0 , Δ t ν 2 I 3 ) and computing

X t j = exp M Δ t A X t j − 1 + V j .

The projected system process p r ( X t ) = X t ∥ X t ∥ can be seen in Figure 8 in blue.

At even time intervals τ 1 , … , τ L = 1 with L = 20 we simulate measurements of X τ j as follows: For each τ j we pseudo-randomly draw W j from N ( 0 , ξ 2 I 3 ) and then compute Z τ j = p r ( pr ( X τ j ) + W j ). These measurements can be seen in Figure 8 in red. The predictions filtered that arise from the application of Algorithm 2 to the measurements Z τ j are visualized in Figure 8 in black.

To illuminate this simulation into something more easily readable, we have also provided Figure 9a which shows the x −coordinate of the same data as in Figure 8 over time but including a 95 %-confidence interval around the filtered mean. Similarly, the y-coordinate is shown in Figure 9b and the z-coordinate in Figure 9c.

As can be seen in Figures 8, 9a–c, the filtered mean is on average closer to the underlying process compared to the distance from the measurements to the underlying process. This tendency will be made clearer in the next set of simulations. We redo the simulation above for ν 2 = 0.1 , 0.2 , 0.5 , 1.0. Then, for each such ν 2 we simulate the filter with 100-fold repetitions for different choices of ξ 2. This is done for

ξ 2 = ν 2 , ν 2 2.8 , ν 2 4.6 , ν 2 6.4 , ν 2 8.2 , ν 2 10 .

In Figures 10a–d, we plot, for different signal to noise ratios in units decibel, i. e 10 ⁡ log 10 ν 2 ξ 2 (dB), the average distance from true point X t to the measurements, and to the filtered means, respectively. As can be seen for all Figures in Figure 10 there is a significant improvement on the average error of the filtered mean compared to the error of the raw measurements. The proposed algorithm (here referred to as the intrinsic EKF) is compared to a naive extrinsic EKF and a second order EKF (in short SOEKF) implemented analogously to Sorenson and Stubberud (1968). All the values of the mean errors in Figure 10 is also tabulated in Table 1.

TABLE 1. Tabulated data corresponding to Figure 10.

5 Simulations on the 4 by 2 Stiefel manifold

Here we consider the filtering problem in Equation 6 on S t 4,2. We choose system parameters ν 2 = 1 and ξ 2 = 0.1, with initial meanand initial covariance Σ 0 = 0.1 i d R 4 × 2.

A = 0 0.173 0.267 − 0.288 − 0.173 0 − 0.279 0.122 − 0.267 0.279 0 0.316 0.288 − 0.122 − 0.316 0 (10)

μ 0 = 1 0 0 1 0 0 0 0 , (11)

The system process X t ∈ R 4,2 is done equivalently to how it is done for S 2 in the previous section. Each coordinate of the projected system process p r ( X t ) = X t ( X t T X t ) − 1 / 2 can be seen in Figure 11 in blue.

Similarly to before, given even time intervals τ 1 , … , τ L = 1 we simulate measurements of X τ j as follows: For each τ j we pseudo-randomly draw W j from N ( 0 . ξ 2 i d R 4 × 2 ) then we compute Z τ j = p r ( pr ( X τ j ) + W j ), here L = 20. These measurements can be seen in red in Figure 11. By directly applying Algorithm 2 with the measurements Z τ j the filtered predictions is visualized in green in Figure 11. Note that here we roughly give 95-% confidence bounds as intervals of the form [ μ K − 1.96 P K , μ K − 1.96 P K ].

As can be seen in Figure 11, on average, the filtered mean is closer to the underlying process compared to the measurements. Again, we shall verify this in the following simulations. We once more define SNR as 10 ⁡ log 10 ν 2 ξ 2 (dB). In Figure 12 one can see the average error of the filtered mean and the average error of the measurements. This is done with the same choices of ν and ξ as was done for S 2 in the previous section. As can be seen in Figures 12a–d, the error of the filtered mean is significantly smaller compared to the error of the measurements. Again, the proposed algorithm (the intrinsic EKF) is compared to a naive extrinsic EKF and a second order EKF (SOEKF). All the values of the mean errors in Figure 12 is also tabulated in Table 2.

TABLE 2. Tabulated data corresponding to Figure 12.

6 Application to a scenario of spacecraft attitude determination

Here we shall give a physically grounded simulation of the orientation of a spacecraft in orbit around the Earth with inspiration following Markley (1998). For example, given the direction of a star, x 1 ∈ R 3 and the (linearly independent) direction of Earth’s magnetic field x 2 ∈ R 3. It is possible to determine the attitude (meaning orientation) of a spacecraft relative to a fixed coordinate system on Earth. For technical reasons, it is only the orthogonal 2 directions spanned by x 1 , x 2 that is necessary in order to determine the attitude of the spacecraft. Therefore, the spacecraft performs measurements that are of the type

p r x 1 x 2 ∈ St 3,2 .

Moreover, as the spacecraft is in orbit of the Earth, the coordinate system is continuously rotated relative to a reference frame on Earth. Therefore, since x 1 and x 2 are fixed points in the fixed reference frame, the observed point p r ( x 1 , X 2 ) ∈ St 3,2 is being continuously rotated by that same rotational motion. This rotational motion is exactly the 1-parameter Lie group generated by the system parameter A in Equation 1. Measurements are performed at the frequency of once every 10 s. If the orbital time for the spacecraft is at about 90 min, and if the unit of time for t is in seconds, we then considerand initial mean

A = 0 0.0012 0 − 0.0012 0 0 0 0 0

μ 0 = 1 0 0 1 0 0 .

We consider uncertainty quantities ν 2 = 1 0 − 8 and ξ 2 = 1 0 − 3, and initial uncertainty Σ 0 = 1 0 − 5 i d R ( 3 × 2 ). Running Algorithm 2 for this scenario over 100 repetitions gives a mean error similar to the extrinsic second order Kalman filter, see Table 3.

TABLE 3. The mean error for the Intrinsic EKF given by Algorithm 2, the extrinsic second order extended Kalman filter (SOEKF) and measurements for the attitude determination of a spacecraft.

7 Discussion

As can be seen in Figures 10 and 12 depending on the situation, one can expect at optimal conditions about a 60% filtering gain. Note, however, that if ν is too large, then the projected distribution will be close to the uniform distribution. In this case the extended Kalman filtering algorithm is not feasible as the predictive power of the dynamics becomes unreliable. This is also true for the error of the measurements, if the error of the measurements are large. Indeed, then there is little information gained from the measurements, naturally.

Moreover, the error of the proposed filtering algorithm seem to be equal in precision to the second order as can be seen in Figures 10 and 12. The proposed algorithm has the upside that no second order derivative needs to be computed, which can be costly to compute as the dimensions go up. Also, for large noise there seem to be little gain of using the proposed algorithm instead of the more simple classical EKF. This is expected, as with large noise the observations does not capture the curvature of the underlying manifold but rather seem like observations from the ambient space (consider, for example S 3 where large noise throws observations around the ball).

Some inaccuracies might originate from the logarithm map used in step 5 of Algorithm 2. The logarithm map is only defined for points close to the base point, therefore one can not guarantee that the algorithm is reliable if the measurements appear outside what is called the injectivity radius (the radius which guarantees that the logarithm map is well defined). It is still an open question what the injectivity radius is for the Stiefel manifold. However, very recently there has been some advances in the following direction: there is a lower bound 4 5 π for the injectivity radius of the Stiefel manifold (under the canonical metric), see Stoye and Zimmermann (2024). Therefore, as long as points are never farther away from each other than 4 5 π, the injectivity of the logarithm map is not an issue.

One obvious limitation with Algorithm 2 is that it requires quite restrictive assumptions on the noise, i.e. we require that the noise is isotropic. Indeed, already on S 2, if we allow for general covariance matrices Σ then the statement of Conjecture 2.4 is rarely true. Thus, if the filtering problem is set up so that future points may not allow for one-to-one correspondence between the covariance of a normal random variable and it’s projected counterpart, our proposed filter is useless. In the case of the sphere, as long as the covariance matrix Σ of the ambient normal random variable has the mean as eigenvector, then Algorithm 2 may still be tractable.

One more subtle issue with Algorithm 2 is that the update steps are only accurate as long as the actual process X t does not significantly drift away from it’s starting point in the radial direction. By choosing A to be anti-symmetric, exp M ( t A ) ∈ S O ( n ), then it is always true that exp M ( t A ) μ 0 ∈ St n , k if μ 0 ∈ S t n , k. Therefore on average, the process X t does not drift away along a radial direction. However, that does not mean that it’s realization does not drift away. In order to bootstrap the filtering process, one way periodically measure the process noise term η ( ν 2 ), of p r ( X t ), since η ( ν 2 ) also depends on the symmetric part S of the polar decomposition X t = p r ( X t ) S.

Finally let us mention some ideas for future improvements. Firstly, observe that the true distribution of p X t ℓ ∣ Z ℓ does not follow the projected normal distribution (though it is close in the L 1-norm) but it is more complicated. Using the true distribution of p X t ℓ ∣ Z ℓ would allow one to adapt the function η to the true distribution that should give performance, although such η might be difficult to attain. Secondly, studying other than Stiefel manifolds would be an interesting direction of future research. For example, one could consider other homogeneous spaces including Grassmann manifolds or more exotic ones such as higher dimensional tori. The latter appears as a model for collection of angles in a protein chain or robotic motion.

Generative AI statement

The author(s) declared that generative AI was used in the creation of this manuscript. During the preparation of this manuscript, the author used Google Gemini (version 3 Pro, January 2026) to assist with literature summarization and text coherence. The author reviewed and edited the output and takes full responsibility for the content.

Any alternative text (alt text) provided alongside figures in this article has been generated by Frontiers with the support of artificial intelligence and reasonable efforts have been made to ensure accuracy, including review by the authors wherever possible. If you identify any issues, please contact us.

Publisher’s note

All claims expressed in this article are solely those of the authors and do not necessarily represent those of their affiliated organizations, or those of the publisher, the editors and the reviewers. Any product that may be evaluated in this article, or claim that may be made by its manufacturer, is not guaranteed or endorsed by the publisher.

Footnotes

Edited by: Congzheng Han, University of Chinese Academy of Sciences, China

Reviewed by: Mahesh Kumar Pal, Indian Institute of Technology Roorkee, India

Junpeng Zhang, Shenzhen Institute of Information Technology, China

Footnotes

This is not an isometric embedding into R n × k when equipped with the Frobenius norm, see Remark 2.3.