See discussions, stats, and author profiles for this publication at: https://www.researchgate.

net/publication/23876857

30 Years of Wahba's Problem

Article · February 1999

Source: NTRS

CITATIONS READS
3 4,553

1 author:

Landis Markley
retired
177 PUBLICATIONS   8,190 CITATIONS   

SEE PROFILE

Some of the authors of this publication are also working on these related projects:

retirement View project

All content following this page was uploaded by Landis Markley on 26 February 2014.

The user has requested enhancement of the downloaded file.

 QUATERNION ATTITUDE ESTIMATION
USING VECTOR OBSERVATIONS
F. Landis Markley* and Daniele Mortari†

ABSTRACT
This paper contains a critical comparison of estimators minimizing Wahba’s
loss function. Some new results are presented for the QUaternion ESTimator
(QUEST) and EStimators of the Optimal Quaternion (ESOQ and ESOQ2) to
avoid the computational burden of sequential rotations in these algorithms. None
of these methods is as robust in principle as Davenport’s q method or the
Singular Value Decomposition (SVD) method, which are significantly slower.
Robustness is only an issue for measurements with widely differing accuracies,
so the fastest estimators, the modified ESOQ and ESOQ2, are well suited to
sensors that track multiple stars with comparable accuracies. More robust forms
of ESOQ and ESOQ2 are developed that are intermediate in speed.

INTRODUCTION
In many spacecraft attitude systems, the attitude observations are naturally represented as unit vectors.
Typical examples are the unit vectors giving the direction to the sun or a star and the unit vector in the
direction of the Earth’s magnetic field. This paper will consider algorithms for estimating spacecraft attitude
from vector measurements taken at a single time, which are known as “single-frame” methods or “point”
methods, instead of filtering methods that employ information about spacecraft dynamics. Almost all
single-frame algorithms are based on a problem proposed in 1965 by Grace Wahba1. Wahba’s problem is to
find the orthogonal matrix A with determinant +1 that minimizes the loss function

∑ab
2
L ( A) ≡ 1
2 i i i − Ari . (1)

where {bi} is a set of unit vectors measured in a spacecraft’s body frame, {ri} are the corresponding unit
vectors in a reference frame, and {ai} are non-negative weights. In this paper we choose the weights to be
inverse variances, ai = σ i−2 , in order to relate Wahba’s problem to Maximum Likelihood Estimation2. This
choice differs from that of Wahba and many other authors, who assumed the weights normalized to unity.
It is possible and has proven very convenient to write the loss function as
L( A) = λ 0 − tr( AB T ) (2)
with
λ 0 ≡ ∑i ai (3)
and
B ≡ ∑i ai bi riT . (4)

Now it is clear that L(A) is minimized when the trace, tr(ABT), is maximized.
* Guidance, Navigation, and Control Center, NASA’s Goddard Space Flight Center, Greenbelt, MD,
landis.markley@gsfc.nasa.gov
† Aerospace Engineering School, University of Rome, “La Sapienza,” Rome, Italy,
daniele@psm.uniroma1.it, Visiting Associate Professor, Texas A&M University, College Station, TX
This has a close relation to the orthogonal Procrustes problem, which is to find the orthogonal matrix A
that is closest to B in the Frobenius norm (also known as the Euclidean, Schur, or Hilbert-Schmidt norm)3
≡ ∑i , j Mij2 = tr( MM T ) .
2
M F
(5)
Now
A− B = A F + B F − 2 tr( AB T ) = 3 + B F − 2 tr( AB T ) ,
2 2 2 2
F
(6)

so Wahba’s problem is equivalent to the orthogonal Procrustes problem with the proviso that the
determinant of A must be +1.
The purpose of this paper is to give an overview in a unified notation of algorithms for solving Wahba’s
problem, to provide accuracy and speed comparisons, and to present two significant enhancements of
existing methods. The popular QUaternion EStimator (QUEST) and EStimators of the Optimal Quaternion
(ESOQ and ESOQ2) algorithms avoid singularities by employing a rotated reference system. Methods
introduced in this paper use information from an a priori quaternion estimate or from the diagonal elements
of the B matrix to determine a desirable reference system, avoiding expensive sequential computations.
Also, tests show that a first-order expansion in the loss function is adequate, avoiding the need for iterative
refinement of the loss function, and motivating the introduction of new first-order versions of ESOQ and
ESOQ2, which are at present the fastest known first-order methods for solving Wahba’s problem.
FIRST SOLUTIONS OF WAHBA’S PROBLEM
J. L. Farrell and J. C. Stuelpnagel 4, R. H. Wessner5, J. R. Velman6, J. E. Brock7, R. Desjardins, and
Wahba presented the first solutions of Wahba’s problem. Farrell and Stuelpnagel noted that any real square
matrix, including B, has the polar decomposition
B = WR, (7)
where W is orthogonal and R is symmetric and positive semidefinite. Then R can be diagonalized by
R = VDV T , (8)
where V is orthogonal and D is diagonal with elements arranged in decreasing order. The optimal attitude
estimate is then given by
Aopt = WV diag[1 1 detW] VT. (9)
In most cases, detW is positive and Aopt = W, but this is not guaranteed.
Wessner and Brock independently proposed the alternative solution
Aopt = (BT)–1(BTB)1/2 = B(BTB)–1/2; (10)
but the matrix inverses in Eq. (10) exist only if B is non-singular, which means that a minimum of three
vectors must be observed. It is well known that two vectors are sufficient to determine the attitude; and the
method of Farrell and Stuelpnagel, as well as the other methods described in this paper, only require B to
have rank two.
SINGULAR VALUE DECOMPOSITION (SVD) METHOD
This method has not been widely used in practice, because of its computational expense, but it yields
valuable analytic insights8,9. The matrix B has the Singular Value Decomposition3:
B = UΣVT = U diag[Σ11 Σ22 Σ33] VT, (11)

where U and V are orthogonal, and the singular values obey the inequalities Σ11 ≥ Σ22 ≥ Σ33 ≥ 0. Then

tr( AB T ) = tr( A V diag[ Σ11 Σ 22 Σ 33 ] U T ) = tr(U T A V diag[ Σ11 Σ 22 Σ 33 ] ) . (12)

The trace is maximized, consistent with the constraint det A = 1, by
U T Aopt V = diag[1 1 (det U )(det V )] , (13)

which gives the optimal attitude matrix

Aopt = U diag[1 1 (det U )(det V )] V T . (14)

The SVD solution is completely equivalent to the original solution by Farrell and Stuelpnagel, since
Eq. (14) is identical to Eq. (9) with U = WV. The difference is that robust SVD algorithms exist now3,10.
In fact, computing the SVD is one of the most robust numerical algorithms.
It is convenient to define
s1 ≡ Σ 11, s2 ≡ Σ 22, and s3 ≡ (detU)(detV) Σ 33, (15)

so that s1 ≥ s2 ≥ |s3|. We will loosely refer to s1, s2, and s3 as the singular values, although the third singular
value of B is actually |s3|. It is clear from Eq. (11) that redefinition of the basis vectors in the reference or
body frame affects V or U, respectively, but does not affect the singular values.
The estimation error is characterized by the rotation angle error vector φerr in the body frame, defined by

exp[φerr ×] = Atrue Aopt

T
, (16a)

where [φ×] is the cross product matrix:

 0 −φ3 φ2 
[φ×] ≡  φ3 0 −φ1  . (16b)
−φ 2 φ1 0 

The SVD method gives the covariance of the rotation angle error vector as
P = U diag[( s2 + s3 ) −1 ( s3 + s1 ) −1 ( s1 + s2 ) −1 ] U T . (16c)

DAVENPORT’S q METHOD
Davenport provided the first useful solution of Wahba’s problem for spacecraft attitude determination11,12.
He parameterized the attitude matrix by a unit quaternion13,14

 q  e sin(φ / 2)
q= = , (17)
q4   cos(φ / 2) 
as
A(q ) = (q42 − q ) I + 2qq T − 2 q4 [q ×] .
2
(18)

The rotation axis e and angle φ will be useful later. Since A(q) is a homogenous quadratic function of q, we
can write
tr( AB T ) = q T Kq (19)
where K is the symmetric traceless matrix
S − I trB z 
K≡
trB
T , (20)
 z
with
 B23 − B32 
S ≡ B+ B T
and z ≡  B31 − B13  = ∑i ai bi × ri . (21)
 B12 − B21 
The optimal attitude is represented by the quaternion maximizing right side of Eq. (19), subject to the unit
constraint q = 1, which is implied by Eq. (17). It is not difficult to see that the optimal quaternion is equal
to the normalized eigenvector of K with the largest eigenvalue, i.e., the solution of
Kqopt ≡ λ max qopt . (22)

With Eqs. (2) and (19), this gives the optimized loss function as
L( Aopt ) = λ o − λ max . (23)

Very robust algorithms exist to solve the symmetric eigenvalue problem3,10.

The eigenvalues of the K matrix, λ max ≡ λ1 ≥ λ 2 ≥ λ3 ≥ λ 4 ≡ λ min , are related to the singular values by11,15
λ 1 = s1 + s2 + s3, λ 2 = s1 – s2 – s3, λ 3 = – s1 + s 2 – s 3 , λ4 = – s1 – s 2 + s 3 . (24)
The eigenvalues sum to zero because K is traceless. There is no unique solution if the two largest
eigenvalues of K are equal, or s2 + s3 = 0. This is not a failure of the q method; it means that the data aren’t
sufficient to determine the attitude uniquely. Equation (16c) shows that the covariance is infinite in this
case. This is expected, since the covariance should be infinite when the attitude is unobservable.
QUATERNION ESTIMATOR (QUEST)
This algorithm, first applied in the MAGSAT mission in 1979, has been the most widely used algorithm
for Wahba’s problem16,17. Equation (22) is equivalent to the two equations
[(λ max + trB) I − S ]q = q4 z (25)
and
(λ max − tr B)q4 = q T z. (26)
Equation (25) gives
q = q4 [(λ max + trB) I − S ]−1 z = {q4 det[(λ max + trB) I − S ]} adj [(λ max + trB) I − S ] z . (27)

The optimal quaternion is then given by

1x 
qopt =  , (28)
γ + x γ 
2 2

where
x ≡ adj [(λ max + tr B) I − S] z = [α I + (λ max − tr B) S + S 2 ] z (29)
and
γ ≡ det[(λ max + trB) I − S ] = α (λ max + tr B) − det S , (30)
with
α ≡ λ2max − ( tr B) 2 + tr (adjS ) . (31)

The second form on the right sides of Eqs. (29) and (30) follows from the Cayley-Hamilton Theorem3,17.
These computations require knowledge of λmax, which is obtained by substituting Eqs. (28) and (29) into
Eq. (26), yielding:
0 = ψ (λ max ) ≡ γ ( λ max − trB) − z T [α I + (λ max − tr B) S + S 2 ] z . (32)

Substituting α and γ from Eqs. (30) and (31) gives a fourth-order equation in λmax, which is simply the
characteristic equation det(λ max I − K ) = 0. Shuster observed that λmax can be easily obtained by Newton-
Raphson iteration of Eq. (32) starting from λ0 as the initial estimate, since Eq. (23) shows that λmax is very
close to λ0 if the optimized loss function is small. In fact, a single iteration is generally sufficient. But
numerical analysts know that solving the characteristic equation is an unreliable way to find eigenvalues, in
general, so QUEST is in principle less robust than Davenport’s original q method. The analytic solution of
the quartic characteristic equation is slower and no more accurate than the iterative solution.
Shuster provided the first estimate of the covariance of the rotation angle error vector in the body frame,

[∑ a ( I − b b )]
−1
P= i i
T
i
i . (33)
He also showed that 2L( Aopt ) obeys a χ probability distribution with 2 nobs − 3 degrees of freedom, to a
2

good approximation and assuming Gaussian measurement errors, where nobs is the number of vector
observations. This can often provide a useful data quality check, as will be seen below.
The optimal quaternion is not defined by Eq. (28) if γ 2 + x = 0 , so it is of interest to see when this
2

condition arises. Applying the Cayley-Hamilton theorem twice to eliminate S4 and S3 after substituting
Eq. (29) gives, with some tedious algebra,
γ 2 + x = γ (dψ / dλ ) ,
2
(34)

where ψ (λ ) is the quartic function defined implicitly by Eq. (32). The discussion following Eq. (15)
implies that dψ / dλ is invariant under rotations, since the coefficients in the polynomial ψ (λ ) depend
only on the singular values of B 15. The Newton-Raphson iteration for λmax requires dψ / dλ to be nonzero,
so γ 2 + x = 0 implies that γ = 0. This means that (qopt ) 4 = 0 and the optimal attitude represents a
2

180˚ rotation. Shuster devised the method of sequential rotations to avoid this singular case16–18.
REFERENCE FRAME ROTATIONS
The (qopt ) 4 = 0 singularity occurs because QUEST does not treat the four components of the quaternion on
an equal basis. Davenport’s q method avoids this singularity by treating the four components
symmetrically, but some other methods have singularities similar to that in QUEST. These singularities
can be avoided by solving for the attitude with respect to a reference coordinate frame related to the original
frame by 180° rotations about the x, y, or z coordinate axis. That is, we solve for one of the quaternions
e i   q  e i   q 4 e i − q × e i 
qi ≡ q ⊗   =   ⊗   =   for i = 1, 2, 3, (35)
 0  q4   0   − q ⋅ e i 
where ei is the unit vector along the ith coordinate axis. We use the convention of Reference 14 for
quaternion multiplication, rather than the historic convention. The products in Eq. (35) are trivial to
implement by merely permuting and changing signs of the quaternion components. For example,
q1 = [q1 , q2 , q3 , q4 ]T ⊗ [1, 0, 0, 0]T = [q4 , − q3 , q2 , − q1 ]T . (36)

The equations for the inverse transformations are the same, since a 180° rotation in the opposite direction
has the same effect. These rotations are also easy to implement on the input data, since a rotation about
axis i simply changes the signs of the jth and kth columns of the B matrix, where {i, j, k} is a permutation
of {1, 2, 3}. The reference system rotation is easily “undone” by Eq. (36) or its equivalent after the optimal
quaternion has been computed.
The original QUEST implementation performed sequential rotations one axis at a time, until an acceptable
reference coordinate system was found. It is clearly preferable to save computations by choosing a single
desirable rotation as early in the computation as possible. This can be accomplished by considering the
components of an a priori quaternion, which is always available in a star tracker application since an a
priori attitude estimate is needed to identify the stars in the tracker’s field of view. If the fourth component
of the a priori quaternion has the largest magnitude, no rotation is performed, while a rotation about the ith
axis is performed if the ith component has the largest magnitude. Then Eq. (36) or its equivalent shows
that the fourth component of the rotated quaternion will have the largest magnitude. This magnitude must
be at least 1/2, but no larger magnitude can be guaranteed, because a unit quaternion may have all four
components with magnitude 1/2. The use of a previous estimate as the a priori attitude guarantees
q4 > cos 82.5° = 0.13 in the rotated frame if rotations between successive estimates are less than 45°.
FAST OPTIMAL ATTITUDE MATRIX (FOAM)
2
The singular value decomposition of B gives a convenient representation for adjB, detB, and B F . These
can be used to write the optimal attitude matrix as15
Aopt = (κλ max − det B) −1[(κ + B F ) B + λ max adjB T − BB T B] ,
2
(37)

where
κ ≡ 12 (λ2max − B F ) .
2
(38)

It’s important to note that all the quantities in Eqs. (37) and (38) can be computed without performing the
singular value decomposition. In this method, λ max is found from

[(
λ max = tr( Aopt B T ) = (κλ max − det B) −1 κ + B
2
F )B 2
F ( )]
+ 3λ max det B − tr BB T BB T , (39)

or, after some matrix algebra,

0 = ψ (λ max ) ≡ (λ2max − B F ) 2 − 8λ max det B − 4 adjB F .
2 2
(40)

Equations (32) and (40) for ψ (λ max ) would be numerically identical with infinite-precision computations,
but the FOAM form of the coefficients is less subject to errors arising in finite-precision computations.
The FOAM algorithm gives the error covariance as:
P = (κλ max − det B) −1 (κ I + BB T ) . (41)

A quaternion can be extracted from Aopt, with a cost of 13 MATLAB flops. This has several advantages: the
four-component quaternion is more economical than the nine-component attitude matrix, easier to
interpolate, and more easily normalized if Aopt is not exactly orthogonal due to computational errors19.
ESTIMATOR OF THE OPTIMAL QUATERNION (ESOQ or ESOQ1)
Davenport’s eigenvalue equation, Eq. (22), says that the optimal quaternion is orthogonal to all the
columns of the matrix
H ≡ K − λ max I , (42)

which means that it must be orthogonal to the three-dimensional subspace spanned by the columns of H.
The optimal quaternion is conveniently computed as the generalized four-dimensional cross-product of any
three columns of this matrix20–22.
Another way of seeing this result is to examine the classical adjoint of H. Representing K in terms of its
eigenvalues and eigenvectors and using the orthonormality of the eigenvectors gives, for any scalar λ,

 4  4
adj( K − λI ) = adj ∑ (λ µ − λ )qµ qµT  = ∑ (λν − λ )(λ ρ − λ )(λτ − λ )qµ qµT . (43)
 µ =1  µ =1

We use Greek indices to label different quaternions, to avoid confusion with Latin indices that label
quaternion components; and let {µ, ν, ρ, τ} denote a permutation of {1, 2, 3, 4}. Setting λ = λ max ≡ λ1
causes all the terms in this sum to vanish except the first, with the result
adj H = (λ 2 − λ max )(λ3 − λ max )(λ 4 − λ max )qopt qopt
T
. (44)

Thus qopt can be computed by normalizing any non-zero column of adj H , which we denote by index k. Let
F denote the symmetric 3×3 matrix obtained by deleting the kth row and kth column from H, and let f
denote the three-component column vector obtained by deleting the kth element from the kth column of H.
Then the kth element of the optimal quaternion is given by
(qopt ) k = − c det F , (45)
and the other three elements are
(qopt )1,L, k −1, k +1,L, 4 = c (adj F )f , (46)

where the coefficient c is determined by normalizing the quaternion. It is desirable to let k denote the
column with the maximum Euclidean norm, which Eq. (44) shows to be the column containing the
maximum diagonal element of the adjoint. Computing all the diagonal elements of adjH, though not as
burdensome as QUEST’s sequential rotations, is somewhat expensive; but this computation can be avoided
by using an a priori quaternion as in QUEST. In the ESOQ case, however, no rotation is performed; we
merely choose k to be the index of the element of the a priori quaternion with maximum magnitude.
The original formulation of ESOQ used the analytic solution of the characteristic equation 23; but the
analytic formula sometimes gives complex eigenvalues, which is theoretically impossible for a real
symmetric matrix. These errors arise from inaccurate values of the coefficients of the quartic characteristic
equation, not from the solution method. It is faster, and equally accurate, to compute λmax by iterative
solution of Eq. (40). Equation (32) would give a faster solution, but it would be less robust, and an even
more efficient solution is described below.
First Order Update (ESOQ1.1)
Test results show that higher-order updates do not improve the performance of the iterative methods,
providing motivation for developing a first-order approximation. The matrix H can be expanded to first order
in ∆λ ≡ λ 0 − λ max as
H = H 0 + ( ∆λ ) I , (47)
where
H 0 ≡ K − λ0 I . (48)

The vector f does not depend on λ max , which only appears in the diagonal elements of H; but the matrix F
depends on λmax, giving

F = F 0 + ( ∆λ ) I , (49)
where F0 is derived from H0 in the same way that F is derived from H. Matrix identities give
adj F = adj F 0 + ∆λ[( tr F 0 ) I − F 0 ] (50)
and
det F = det F 0 + ( ∆λ ) tr(adj F 0 ) , (51)

to first order in ∆λ . The characteristic equation can be expressed to the same order as

[ ]
0 = det H = ( Hkk0 + ∆λ ) det F − f T (adj F ) f = Hkk0 detF 0 − f T g + Hkk0 tr(adj F 0 ) + det F 0 − f T h ∆λ , (52)

where
g ≡ (adjF 0 )f . (53)

and
h ≡ [( trF 0 ) I − F 0 ]f . (54)
Equation (52) is easily solved for ∆λ ≡ λ 0 − λ max , and then the first order quaternion estimate is given by

(qopt ) k = − c[det F 0 + ( ∆λ ) tr(adj F 0 )] (55)

and
(qopt )1,L, k −1, k +1,L, 4 = c (g + h ∆λ ) . (56)
SECOND ESTIMATOR OF THE OPTIMAL QUATERNION (ESOQ2)
This algorithm uses the rotation axis/angle form of the optimal quaternion, as given in Eq. (17).
Substituting these into Eqs. (25) and (26) gives
(λ max − tr B) cos(φ / 2) = z T e sin(φ / 2) (57)
and
z cos(φ / 2) = [(λ max + trB) I − S ] e sin(φ / 2) (58)

Multiplying Eq. (58) by (λ max − tr B) and substituting Eq. (57) gives

Me sin(φ / 2) = 0 , (59)
where
M ≡ (λ max − trB)[(λ max + trB) I − S ] − zz T = [m1 M m 2 M m 3 ] . (60)

These computations lose numerical significance if (λ max − tr B) and z are both close to zero, which would
be the case for zero rotation angle. We can always avoid this singular condition by using one of the
sequential reference system rotations16–18 to ensure that trB is less than or equal to zero. If we rotate the
reference frame about the ith axis
( tr B) rotated = ( Bii − Bjj − Bkk ) unrotated = (2 Bii − tr B) unrotated , (61)
where {i, j, k} is a permutation of {1, 2, 3}. Thus no rotation is performed in ESOQ2 if trB is the
minimum of {B11 , B22 , B33 , tr B} , while a rotation about the ith axis is performed if Bii is the minimum.
This will ensure the most negative value for the trace in the rotated frame. As in the QUEST case, the
rotation is easily “undone” by Eq. (36) or its equivalent after the quaternion has been computed. Note that
efficiently finding an acceptable rotated frame for ESOQ2 does not require an a priori attitude estimate.
Equation (59) says that the rotation axis is a null vector of M. The columns of adj M are the cross products
of the columns of M:
adj M = [m 2 × m 3 M m 3 × m1 M m1 × m 2 ] . (62)

Because M is singular, all these columns are parallel, and all are parallel to the rotation axis e. Thus we set
e=y y, (63)

where y is the column of adj M (i.e., the cross product) with maximum norm. Because M is symmetric, it
is only necessary to find the maximum diagonal element of its adjoint to determine which column to use.
The rotation angle is found from Eq. (57) or one of the components of Eq. (58). We will show that Eq. (57)
is the best choice. Comparing Eq. (20) with the eigenvector/eigenvalue expansion
h
K = ∑ λ µ qµ qµT , (64)
µ =1

establishes the identities

4
tr B = ∑ λ µ (qµ ) 24 , (65)
µ =1

and
4
z = ∑ λ µ (qµ ) 4 q µ . (66)
µ =1

Using Eq. (65) and the orthonormality of the eigenvectors of K, we find that
4
z = ∑ λ2µ (qµ ) 24 − (tr B) 2 .
2
(67)
µ =1
 4
This equation and ∑ (qµ ) 2
4 = 1 give the inequality
µ =1

z ≤ max λ µ = max (λ max , − λ min ) . (68)

µ =1,L,4

This shows that choosing the rotated reference system that provides the most negative value of trB makes
Eq. (57) the best equation to solve for the rotation angle. With Eq. (63), this can be written
(λ max − tr B) y cos(φ / 2) = (z ⋅ y)sin(φ / 2) , (69)

which means that there is some scalar η for which

cos(φ / 2) = η (z ⋅ y) (70)

and
sin(φ / 2) = η (λ max − tr B) y . (71)

Substituting into Eq. (17) and using Eq. (63) gives the optimal quaternion as24,25
1 (λ max − tr B) y 
qopt =  . (72)
(λ max − tr B) y + (z ⋅ y) 
2 2 z⋅y 

Note that it is not necessary to normalize the rotation axis. ESOQ2 does not define the rotation axis
uniquely if M has rank less than two. This includes the usual case of unobservable attitude and also the case
of zero rotation angle. Requiring trB to be non-positive avoids zero rotation angle singularity, however. We
compute λmax by iterative solution of Eq. (40) in the general case, as for ESOQ.

First Order Update (ESOQ2.1)

The motivation for and development of this algorithm are similar to those of the first order update used in
ESOQ1.1. The matrix M can be expanded to first order in ∆λ ≡ λ 0 − λ max as

M = M 0 + ( ∆λ ) N , (73)
where
M 0 ≡ (λ 0 − trB)[(λ 0 + trB) I − S ] − zz T = [m10 M m 20 M m 30 ] (74)

and
N ≡ S − 2 λ 0 I = [ n1 M n 2 M n 3 ] . (75)
To this same order, we have
y ≡ m i × m j = (m i0 + n i ∆λ ) × (m 0j + n j ∆λ ) = y 0 + p ∆λ , (76)

where
y 0 ≡ m i0 × m 0j . (77)

and
p ≡ m i0 × n j + n i × m 0j . (78)

The maximum eigenvalue can be found from the condition that M is singular. This gives the first-order
approximation
0 = det M = (m i × m j ) ⋅ m k = (y 0 + p ∆λ ) ⋅ (m 0k + n k ∆λ ) = y 0 ⋅ m 0k + (y 0 ⋅ n k + m 0k ⋅ p)∆λ , (79)

where {i, j, k} is a cyclic permutation of {1, 2, 3}. This is solved for ∆λ , and the attitude estimate is found
by substituting Eq. (76) and λ max = λ 0 − ∆λ into Eq. (72).
There is an interesting relation between the eigenvalue condition det M (λ ) = 0 used in ESOQ2.1 and the
condition ψ (λ ) = 0 used in other algorithms. Since M(λ ) is a 3×3 matrix quadratic in λ, the eigenvalue
condition is of sixth order in λ. Straightforward matrix algebra shows that
det M (λ ) = (λ − tr B) 2 ψ (λ ) . (80)

Thus det M (λ ) has the four roots of ψ (λ ), the eigenvalues of Davenport’s K matrix, and an additional
double root at trB. Choosing the rotated reference frame maximizing – trB ensures that this double root is
far from the desired root at λmax.
TESTS
We test the accuracy and speed of MATLAB implementations of these methods, using simulated data. The q
and SVD methods use the functions eig and svd, respectively; the others use the equations in this paper.
MATLAB uses IEEE double-precision floating-point arithmetic, in which the numbers have approximately
16 significant decimal digits26.
We analyze three test scenarios. In all these scenarios, the pointing of one spacecraft axis, which we take to
be the spacecraft x axis, is much better determined that the rotation about this axis. This is a very common
case that arises in spacecraft that point a single instrument (like an astronomical telescope) very precisely.
This is also a characteristic of attitude estimates from a single narrow-field-of-view star tracker, where the
rotation about the tracker boresight is much less well determined than the pointing of the boresight. The x
axis error and the yz error, which is the error about an axis orthogonal to the x axis and determines the x
axis pointing, are computed from an error quaternion qerr by writing
 q err  e x sin(φ x / 2) e yz sin(φ yz / 2)
qerr =  = ⊗ 
qerr4   cos(φ x / 2)   cos(φ yz / 2) 
(81)
e x cos(φ yz / 2)sin(φ x / 2) + e yz cos(φ x / 2)sin(φ yz / 2) − (e x × e yz )sin(φ x / 2)sin(φ yz / 2)
= ,
 cos(φ x / 2) cos(φ yz / 2) 
where ex = [1 0 0]T and eyz is a unit vector orthogonal to ex. We can always find φx in [–π, π] and φyz in
[0, π] by selecting eyz appropriately. Then, since eyz and ex × eyz form an orthonormal basis for the y-z
(or 2-3) plane, the error angles are given by
φ x = 2 tan −1 (qerr1 qerr4 ) (82)
and
φ yz = 2 sin −1 ( 2
qerr2 + qerr3
2
). (83)

Equations (82) and (83) would be unchanged if the order of the rotations about ex and eyz were reversed; only
the unit vector eyz would be different.
The magnitude φ err of the rotation angle error vector defined by Eq. (16a) is given by

cos(φ err / 2) = qerr4 = cos(φ x / 2) cos(φ yz / 2) . (84)

Thus φ err / 2 is the hypotenuse of a right spherical triangle with sides φ x / 2 and φ yz / 2 . The angles φx and
φyz are the spherical trigonometry equivalents of two orthogonal components of the error vector.
First Test Scenario
The first scenario simulates an application for which the QUEST algorithm has been widely used. A single
star tracker with a narrow field of view and boresight at [1, 0, 0]T is assumed to track five stars at
1  0.99712   0.99712  0.99712   0.99712 
b1 = 0 , b 2 = 0.07584 , b3 = − 0.07584 , b 4 =  0 , and b 5 = 
        0 .
 (85)
0   0   0  0.07584  − 0.07584 
We simulate 1000 test cases with uniformly distributed random attitude matrices, which we use to map the
five observation vectors to the reference frame. The reference vectors are corrupted by Gaussian random
noise with equal standard deviations of 6 arcseconds per axis and then normalized. Equation (33) gives the
covariance for the star tracker scenario as
5
P = (6 arcsec) 2 [5 I − ∑ bi biT ]−1 = diag[1565, 7.2, 7.2] arcsec 2 , (86)
i =1

which gives the standard deviations of the attitude estimation errors as

σ x = 1565 arcsec = 40 arcsec and σ yz = 7.2 + 7.2 arcsec = 3.8 arcsec . (87)

This is a very favorable five-star case, since the stars are uniformly and symmetrically distributed across the
tracker’s field of view. One advantage of simulating a fixed star distribution and applying Gaussian noise to
the reference vectors rather than the observation vectors is that Eqs. (85–87) are always valid, and the
predicted covariance can be compared with the results of the Monte Carlo simulation.

The loss function is computed with measurement variances in (radians)2, since this results in 2L( Aopt )
approximately obeying a χ2 distribution. The minimum and maximum values of the loss function in the
1000 test runs are 0.23 and 12.1, respectively. The probability distribution of the loss function is plotted as
the solid line in Figure 1, and several values of P( χ 2 | ν ) for χ 2 = 2 L( Aopt ) and ν = 2 nobs − 3 = 7 are
plotted as circles23. The agreement is seen to be excellent, which indicates that the measurement weights
accurately reflect the normally-distributed measurement errors for this scenario.
The RSS (outside of parentheses) and the maximum (in parentheses) estimation errors over the 1000 cases
for the star tracker scenario are presented in Table 1. The q method and the SVD method should both give
the truly optimal solution, since they are based on robust matrix analysis algorithms3,10. The q method is
taken as optimal by definition, so no estimated-to-optimal differences are presented for that algorithm, and
the differences between the SVD and q methods provide an estimate of the computational errors of both
methods. No estimate of the loss function is provided when no update of λmax is performed, accounting for
the lack of entries in the loss function column in the tables for these cases.
The loss function is computed exactly by both the q and SVD methods, in principle. Equation (3) gives
λ 0 = 5.9 × 10 9 rad–2 for this scenario, so the expected errors in double-precision machine computation of
λmax and thus of the loss function is on the order of 10 −16 times this, or about 10 −6 , in rough agreement
with the difference shown in the table. In fact, all the algorithms that compute the loss function give
nearly the same accuracy in this scenario. The table also shows show that one Newton-Raphson iteration
for λmax is always sufficient, with a second iteration providing no practical improvement.
Although the attitude estimates of the different algorithms in Table 1 may be closer or farther from the
optimal estimate, all the algorithms provide estimates that are equally close to the true attitude. The number
of digits presented in the table is chosen to emphasize this point, but these digits are not all significant; the
results of 1000 different random cases would not agree with these cases to more than two decimal places.
The smallest angle differences in the table, about 10 −10 arcseconds, are equal to 5 × 10 −16 radians, which is
at the limit of double precision math. The differences between the estimated and optimal values further
show that no update of λmax is required in this scenario, since the estimates using λ0 are equally close to the
truth. Finally, it is apparent that the covariance estimate of Eq. (86) is quite accurate.
 1

0.9

0.8

0.7
probability distribution function

0.6

0.5

0.4

0.3

0.2

0.1

0
0 2 4 6 8 10 12 14
loss function

Figure 1: Empirical (solid line) and Theoretical (circles) Loss Function Distribution
for the Seven-Degree-of-Freedom Star Tracker Scenario
Second Test Scenario
The second scenario uses three observations with widely varying accuracies to provide a difficult test case
for the algorithms under consideration. The three observation vectors are
1  −0.99712   − 0.99712 
 , = 
b1 = 0  b 2  0 07584 
.  , and b3  − 0.07584 
= . (88)
0   0   0 

We simulate 1000 test cases as in the star tracker scenario, but with Gaussian noise of one arcsecond per
axis on the first observation, and 1° per axis on the other two. This models the case that the first
observation is from an onboard astronomical telescope, and the other two observations are from a coarse sun
sensor and a magnetometer, for example. A very accurate estimate of the orientation of the x axis is required
in such an application, but the rotation about this axis is expected to be fairly poorly determined. This is
reflected by the predicted covariance in this scenario, which is to a very good approximation,

P = diag[ 12 (1 − 0.99712 2 ) −1 deg 2 , 1 arcsec 2 , 1 arcsec 2 ] , (89)

giving
σ x = 9.3 deg and σ yz = 1.4 arcsec . (90)
 Table 1: Estimation Errors for Star Tracker Scenario
Algorithm RSS (max) estimated–to–optimal RSS (max) estimated–to–true
λmax iterations loss function x (arcsec) yz (arcsec) x (arcsec) yz (arcsec)
q — — — — 38.41 (122.5) 3.829 (9.252)
SVD — 0.4 (1.8) × 10–5 1.4 (5.6) × 10–8 0.8 (2.9) × 10–10 38.41 (122.5) 3.829 (9.252)
0 — 0.014 (0.078) 9.7 (36) × 10 –5
38.41 (122.5) 3.829 (9.252)
FOAM 1 0.4 (1.6) × 10 –5
1.5 (5.6) × 10 –8
26 (104) × 10 –10
38.41 (122.5) 3.829 (9.252)
2 0.4 (1.7) × 10 –5
1.5 (5.6) × 10 –8
26 (88) × 10 –10
38.41 (122.5) 3.829 (9.252)
0 — 0.015 (0.078) 9.6 (36) × 10–5 38.41 (122.5) 3.829 (9.252)
QUEST 1 2.5 (7.2) × 10 –5
10.1 (46) × 10 –8
6.1 (26) × 10 –10
38.41 (122.5) 3.829 (9.252)
2 2.9 (8.4) × 10 –5
11.1 (50) × 10 –8
7.2 (25) × 10 –10
38.41 (122.5) 3.829 (9.252)
0 — 0.014 (0.078) 9.7 (36) × 10 –5
38.41 (122.5) 3.829 (9.252)
ESOQ2 1 0.4 (1.7) × 10–5 1.5 (6.1) × 10–8 2.0 (10) × 10–10 38.41 (122.5) 3.829 (9.252)
2 0.4 (1.8) × 10 –5
1.5 (6.1) × 10 –8
2.1 (11) × 10 –10
38.41 (122.5) 3.829 (9.252)
ESOQ2.1 1 0.4 (1.6) × 10 –5
1.5 (5.9) × 10 –8
1.9 (12) × 10 –10
38.41 (122.5) 3.829 (9.252)
0 — 0.015 (0.078) 9.6 (36) × 10 –5
38.41 (122.5) 3.829 (9.252)
ESOQ 1 0.4 (1.7) × 10–5 1.5 (6.2) × 10–8 9.6 (39) × 10–10 38.41 (122.5) 3.829 (9.252)
2 0.4 (1.8) × 10 –5
1.5 (6.2) × 10 –8
9.6 (39) × 10 –10
38.41 (122.5) 3.829 (9.252)
ESOQ1.1 1 1.0 (5.3) × 10 –5
4.1 (24) × 10 –8
7.0 (29) × 10 –10
38.41 (122.5) 3.829 (9.252)

The minimum and maximum values of the loss function computed by the q method in the 1000 test runs
for the second scenario are 0.003 and 8.5, respectively. The probability distribution of the loss function is
plotted as the solid line in Figure 2, and several values of the χ 2 distribution with three degrees of freedom
are plotted as circles. The agreement is almost as good as the seven-degree-of-freedom case.
The estimation errors for this scenario are presented in Table 2, which is similar to Table 1 except that the
rotation errors about the x axis are given in degrees. The agreement of the q and SVD methods is virtually
identical to their agreement in the star tracker scenario, but the other algorithms show varying performance.
The best results for the attitude accuracies are in agreement with the covariance estimates of Eq. (89).
Equation (3) gives λ 0 = 8.5 × 1010 rad–2 for this scenario, so the expected accuracy of the loss function in
double-precision machine computations is on the order of 10–5, which is the level of agreement between the
values computed by the q and SVD methods. None of the other methods computes the loss function nearly
as accurately. This differs from the first scenario, where all the algorithms came close to achieving the
maximum precision available in double-precision arithmetic. The iterative computation of λmax in QUEST,
ESOQ1.1, and ESOQ2.1 is poor, but this has surprisingly little effect on the determination of the x axis
pointing. The determination of the rotation about the x axis is adversely affected by an inaccurate
computation of λmax, however, with maximum deviations from the optimal estimate of almost 180°. The
only useful results of QUEST are obtained by not performing any iterations for λmax.
As noted in the discussion of the analytic solution of the characteristic equation for ESOQ1, errors in the
computation of the eigenvalues are believed to arise from inaccurate values of the coefficients of the quartic
characteristic equation rather than from the solution method employed. The superior accuracy of the iterative
computation of λmax in FOAM, ESOQ, and ESOQ2 as compared to QUEST, ESOQ1.1, and ESOQ2.1 is
likely due to the fact that Eq. (40) deals with B directly, while the other algorithms lose some numerical
significance by using the symmetric and skew parts S and z.
 1

0.9

0.8

0.7
probability distribution function

0.6

0.5

0.4

0.3

0.2

0.1

0
0 1 2 3 4 5 6 7 8 9
loss function

Figure 2: Empirical (solid line) and Theoretical (circles) Loss Function Distribution
for the Three-Degree-of-Freedom Unequal Measurement Weight Scenario
Third Test Scenario
The third scenario investigates the effect of measurement noise mismodeling, illustrating problems that first
appeared in analyzing data from the Upper Atmosphere Research Satellite27. Of course, no one would
intentionally use erroneous models, but it can be very difficult to determine an accurate noise model for real
data, and the assumption of any level of white noise is often a poor approximation to real measurement
errors. This scenario uses the same three observation vectors as the second scenario, given by Eq. (88). We
again simulate 1000 test cases, but with Gaussian white noise of 1° per axis on the first observation and
0.1° per axis on the other two. The estimator, however, incorrectly assumes measurement errors of 0.1° per
axis on all three observations, so it weights the measurements equally.
The minimum and maximum values of the loss function computed by the q method in the 1000 test runs
for the third scenario are 0.07 and 453, respectively. The probability distribution of the loss function is
plotted in Figure 3. The theoretical three-degree-of-freedom distribution was plotted in Figure 2; it is not
plotted in Figure 3 since it would be a very poor fit to the data. More than 95% of the values of L( Aopt ) are
theoretically expected to lie below 4, according to the χ 2 distribution of Figure 2, but almost half of the
values of the loss function in Figure 3 have values greater than 50. Shuster has emphasized that large
values of the loss function are an excellent indication of measurement mismodeling or simply of bad data.
 Table 2: Estimation Errors for Unequal Measurement Weight Scenario
Algorithm RSS (max) estimated–to–optimal RSS (max) estimated–to–true
λmax iterations loss function x (deg) yz (arcsec) x (deg) yz (arcsec)
q — — — — 9.5 (34) 1.42 (3.57)
SVD — 1.6 (6.9) × 10 –5
1.4 (8.0) × 10 –5
7.7 (24) × 10 –11
9.5 (34) 1.42 (3.57)
0 — 1.5 (9.9) 7.9 (29) × 10 –3
9.5 (34) 1.42 (3.57)
FOAM 1 0.09 (0.7) 0.09 (1.1) 8.0 (26) × 10–3 9.5 (34) 1.42 (3.57)
2 0.0007 (0.012) 0.0008 (0.013) 7.8 (29) × 10 –3
9.5 (34) 1.42 (3.57)
0 — 1.9 (12) 0.4 (3.8) × 10 –3
9.6 (34) 1.42 (3.57)
QUEST 1 768 (2329) 60 (170) 2.3 (9.0) × 10 –3
48 (90) 1.42 (3.57)
2 1796 (38501) 62 (175) 4.8 (95) × 10–3 48 (91) 1.42 (3.57)
0 — 1.5 (9.9) 1.1 (6.8) × 10 –3
9.5 (34) 1.42 (3.57)
ESOQ2 1 0.09 (0.7) 0.09 (1.1) 1.3 (9.4) × 10 –3
9.5 (34) 1.42 (3.57)
2 0.0007 (0.012) 0.0008 (0.013) 1.1 (7.1) × 10 –3
9.5 (34) 1.42 (3.57)
ESOQ2.1 1 59 (370) 39 (178) 1.5 (14) × 10–3 29 (91) 1.42 (3.57)
0 — 1.9 (12) 4.8 (23) × 10 –3
9.6 (34) 1.42 (3.57)
ESOQ 1 0.09 (0.7) 0.10 (1.1) 5.3 (28) × 10 –3
9.5 (34) 1.42 (3.57)
2 0.0007 (0.012) 0.0008 (0.013) 5.2 (24) × 10 –3
9.5 (34) 1.42 (3.57)
ESOQ1.1 1 327 (1727) 60 (177) 2.6 (34) × 10–3 43 (90) 1.42 (3.57)

The estimation errors for this scenario are presented in Table 3, which is similar to Tables 1 and 2 except
that all the angular errors are given in degrees. The truly optimal q and SVD methods agree even more
closely than in the other scenarios. Equation (3) gives λ 0 = 5 × 10 5 rad–2 for this scenario, so the expected
accuracy of the loss function in double-precision machine computations is on the order of 10–10, the level of
agreement between the q and SVD methods. As in the second scenario, none of the other methods computes
the loss function nearly as accurately. In the third scenario, though, the iterative computation of λmax works
well for all the algorithms, and both iterations improve the agreement of the loss function and attitude
estimates with the optimal values. The first order refinement is reflected in a reduction of the attitude errors,
particularly in determining the rotation about the x axis, but no algorithm is aided significantly by a
second-order update. As in the first scenario, all the algorithms with the first order update to λmax perform as
well as the q and SVD methods.
Speed
Absolute speed numbers are not very important for ground computations, since the actual estimation
algorithm is only a part of the overall attitude determination data processing effort. Speed was more
important in the past, when thousands of attitude solutions had to be computed by slower machines, which
is why QUEST was so important for the MAGSAT mission. For a real-time computer in a spacecraft
attitude control system or a star tracker, which must finish all its required tasks in a limited time, the
longest computation time is more important than the average time. This would penalize some algorithms
for real-time applications, unless we eliminate the need for sequential rotations by using the methods
described above. Our speed comparisons use an a priori quaternion for QUEST and ESOQ, or the diagonal
elements of B for ESOQ2, to eliminate these extra computations. Its independence from a priori attitude
information somewhat favors ESOQ2 for real-time applications.
 1

0.9

0.8

0.7
probability distribution function

0.6

0.5

0.4

0.3

0.2

0.1

0
0 50 100 150 200 250 300 350 400 450 500
loss function

Figure 3: Empirical Loss Function Distribution for the

Mismodeled Measurement Weight Scenario
Figures 4 and 5 show the maximum number of MATLAB floating-point operations (flops) to compute an
attitude using two to six reference vectors; the times to process more than six vectors follow the trends seen
in the figure. The inputs for the timing tests are the nobs normalized reference and observation vector pairs
and their nobs weights. One thousand test cases with random attitudes and random reference vectors with
Gaussian measurement noise were simulated for each number of reference vectors.
Figure 4 plots the times of the more robust methods. The break in the plots for FOAM, ESOQ, and
ESOQ2 at nobs = 3 results from using an exact solution of the characteristic equation in the two-observation
case, when det B = 0 and Eq. (40) shows that ψ (λ max ) is a quadratic function of λ2max . For three or more
observations, these algorithms are timed for a first-order update to λmax using Eq. (40). Additional iterations
for λmax are not expensive, costing only 11 flops each. It is clear that the q method and the SVD method
require significantly more computational effort than the other algorithms, as expected. The q method is
more efficient than the SVD method, except in the least interesting two-observation case. The other three
algorithms are much faster, with the fastest, ESOQ and ESOQ2, being nearly equal in speed. An
implementation of the q method computing only the largest eigenvalue of the K matrix and its eigenvector
would be faster in principle than eig, which computes all four eigenvalues and eigenvectors. However, the
MATLAB routine that can be configured this way is much slower than eig. This option was not
investigated further, since the q method is unlikely to be the method of choice when speed is a primary
consideration.
 Table 3: Estimation Errors for Mismodeled Measurement Weight Scenario
Algorithm RSS (max) estimated–to–optimal RSS (max) estimated–to–true
λmax iterations loss function x (deg) yz (deg) x (deg) yz (deg)
q — — — — 0.96 (3.62) 0.49 (1.17)
SVD — 4.1 (22) × 10–10 3.8 (17) × 10–12 2.3 (7.3) × 10–14 0.96 (3.62) 0.49 (1.17)
0 — 0.7 (5.9) 4.0 (21) × 10 –3
1.18 (5.42) 0.49 (1.16)
FOAM 1 2.6 (24) 0.020 (0.33) 1.0 (11) × 10 –4
0.96 (3.60) 0.49 (1.17)
2 0.004 (0.07) 0.4 (10) × 10 –4
1.7 (35) × 10 –7
0.96 (3.62) 0.49 (1.17)
0 — 0.9 (7.6) 4.3 (29) × 10–3 1.27 (7.66) 0.49 (1.16)
QUEST 1 2.6 (24) 0.023 (0.33) 1.1 (11) × 10 –4
0.96 (3.60) 0.49 (1.17)
2 0.004 (0.07) 0.4 (10) × 10 –4
1.7 (35) × 10 –7
0.96 (3.62) 0.49 (1.17)
0 — 0.7 (5.9) 4.0 (21) × 10 –3
1.18 (5.42) 0.49 (1.16)
ESOQ2 1 2.6 (24) 0.020 (0.33) 1.0 (11) × 10–4 0.96 (3.60) 0.49 (1.17)
2 0.004 (0.07) 0.4 (10) × 10 –4
1.7 (35) × 10 –7
0.96 (3.62) 0.49 (1.17)
ESOQ2.1 1 2.6 (24) 0.020 (0.33) 0.6 (5.8) × 10 –4
0.96 (3.60) 0.49 (1.17)
0 — 0.9 (7.6) 4.3 (29) × 10 –3
1.27 (7.66) 0.49 (1.16)
ESOQ 1 2.6 (24) 0.023 (0.33) 1.1 (11) × 10–4 0.96 (3.60) 0.49 (1.17)
2 0.004 (0.07) 0.4 (10) × 10 –4
1.7 (35) × 10 –7
0.96 (3.62) 0.49 (1.17)
ESOQ1.1 1 2.6 (24) 0.023 (0.33) 1.3 (27) × 10 –6
0.96 (3.60) 0.49 (1.17)

Figure 5 compares the timing of the fastest methods, which generally use zeroth order and first order
approximations for λ max . Both QUEST(1) and ESOQ2.1 use the exact quadratic solution for λ max in the
two-observation case, but ESOQ1.1 uses its faster first order approximation for any number of
observations. It is clear that ESOQ and ESOQ2 are the fastest algorithms using the zeroth order
approximation for λmax, and ESOQ1.1 is the fastest of the first order methods.
CONCLUSIONS
This paper has examined the most useful algorithms for estimating spacecraft attitude from vector
measurements based on minimizing Wahba’s loss function. These were tested in three scenarios, which
show that the most robust, reliable, and accurate estimators are Davenport’s q method and the Singular
Value Decomposition (SVD) method. This is not surprising, since these methods are based on robust and
well-tested general-purpose matrix algorithms. The q method, which computes the optimal quaternion as the
eigenvector of a symmetric 4×4 matrix with the largest eigenvalue, is the faster of these two.
Several algorithms are significantly less burdensome computationally than the q and SVD methods. These
methods are less robust in principle, since they solve the quartic characteristic polynomial equation for the
maximum eigenvalue, a procedure that is potentially numerically unreliable. Algorithms that use the form
of the characteristic polynomial from the Fast Optimal Attitude Matrix (FOAM) algorithm performed as
well as the q and SVD methods in practice, however. The fastest of these algorithms are the EStimators of
the Optimal Quaternion, ESOQ and ESOQ2. The execution times of these methods are reduced by using the
information from an a priori attitude estimate to eliminate sequential rotations in QUEST and extra
computations in ESOQ, or information derived from the observations to speed ESOQ2.
 1200

1100 SVD

1000

q method
900
floating point operations

800

700

600

500

400
FOAM
300
ESOQ2
ESOQ
200
2 3 4 5 6
number of observed vectors

Figure 4: Execution Times for Robust Estimation Algorithms

FOAM, ESOQ, and ESOQ2 use first order approximation for λmax.
All the algorithms tested perform as well as the more robust algorithms in cases where measurement
weights do not vary too widely and are reasonably well modeled. If the measurement uncertainties are not
well represented by white noise, however, an update is required, while this update can be unreliable if the
measurement weights span a wide range. The examples in the paper show that these robustness concerns are
not an issue for the processing of multiple star observations with comparable accuracies, the most common
application of Wahba’s loss function. Thus the fastest algorithms, the zeroth-order ESOQ and ESOQ2 and
the first-order ESOQ1.1, are well suited to star tracker attitude determination applications. In general-
purpose applications where measurement weights may vary greatly, one of the more robust algorithms may
be preferred.
 280

260

240

QUEST(1)
220
floating point operations

ESOQ2.1 ESOQ1.1
200

QUEST(0)
180

ESOQ2(0)
160 ESOQ(0)

140

120

100
2 3 4 5 6
number of observed vectors

Figure 5: Execution Times for Fast Estimation Algorithms

Numbers in parentheses denote order of λmax approximation.

REFERENCES
1. Wahba, Grace, “A Least Squares Estimate of Spacecraft Attitude,” SIAM Review, Vol. 7, No. 3,
July 1965, p. 409.
2. Shuster, Malcolm D., “Maximum Likelihood Estimate of Spacecraft Attitude,”
Journal of the Astronautical Sciences, Vol. 37, No. 1, January-March 1989, pp. 79–88.
3. Horn, Roger A. and Charles R. Johnson, Matrix Analysis, Cambridge, UK,
Cambridge University Press, 1985.
4. Farrell, J. L. and J. C. Stuelpnagel, “A Least Squares Estimate of Spacecraft Attitude,”
SIAM Review, Vol. 8, No. 3, July 1966, pp. 384-386.
5. Wessner, R. H., SIAM Review, Vol. 8, No. 3, July 1966, p. 386.
6. Velman, J. R., SIAM Review, Vol. 8, No. 3, July 1966, p. 386.
7. Brock, J. E., SIAM Review, Vol. 8, No. 3, July 1966, p. 386.
8. Markley, F. Landis, “Attitude Determination Using Vector Observations and the Singular Value
Decomposition,” AAS Paper 87-490, AAS/AIAA Astrodynamics Specialist Conference, Kalispell, MT,
August 1987.
 9. Markley, F. Landis, “Attitude Determination Using Vector Observations and the Singular Value
Decomposition,” Journal of the Astronautical Sciences, Vol. 36, No. 3, July-Sept. 1988, pp. 245-258.
10. Golub, Gene H. and Charles F. Van Loan, Matrix Computations, Baltimore, MD,
The Johns Hopkins University Press, 1983.
11. Keat, J., “Analysis of Least-Squares Attitude Determination Routine DOAOP,”
Computer Sciences Corporation Report CSC/TM-77/6034, February 1977.
12. Lerner, Gerald M., “Three-Axis Attitude Determination,” in Spacecraft Attitude Determination and
Control, ed. by James R. Wertz, Dordrecht, Holland, D. Reidel, 1978.
13. Markley, F. Landis, “Parameterizations of the Attitude,” in Spacecraft Attitude Determination and
Control, ed. by James R. Wertz, Dordrecht, Holland, D. Reidel, 1978.
14. Shuster, Malcolm D., “A Survey of Attitude Representations,” Journal of the Astronautical Sciences,
Vol. 41, No. 4, October-December 1993, pp. 439-517.
15. Markley, F. Landis, “Attitude Determination Using Vector Observations: a Fast Optimal Matrix
Algorithm,” Journal of the Astronautical Sciences, Vol. 41, No. 2, April-June 1993, pp. 261-280.
16. Shuster, M. D. “Approximate Algorithms for Fast Optimal Attitude Computation,”
AIAA Paper 78-1249, AIAA Guidance and Control Conference, Palo Alto, CA, August 7–9, 1978.
17. Shuster, M. D. and S. D. Oh, "Three-Axis Attitude Determination from Vector Observations,"
Journal of Guidance and Control, Vol. 4, No. 1, January-February 1981, pp. 70-77.
18. Shuster, Malcolm D. and Gregory A. Natanson, “Quaternion Computation from a Geometric Point of
View,” Journal of the Astronautical Sciences, Vol. 41, No. 4, October-December 1993, pp. 545-556.
19. Markley, F. L., “New Quaternion Attitude Estimation Method,” Journal of Guidance, Control, and
Dynamics, Vol. 17, No. 2, March-April 1994, pp. 407-409.
20. Mortari, Daniele, “EULER-q Algorithm for Attitude Determination from Vector Observations,”
Journal of Guidance, Control, and Dynamics, Vol. 21, No. 2, March-April 1998, pp. 328-334.
21. Mortari, Daniele, “ESOQ: A Closed-Form Solution to the Wahba Problem,”
Journal of the Astronautical Sciences, Vol. 45, No.2 April-June 1997, pp. 195-204.
22. Mortari, Daniele, “n-Dimensional Cross Product and its Application to Matrix Eigenanalysis,”
Journal of Guidance, Control, and Dynamics, Vol. 20, No. 3, May-June 1997, pp. 509-515.
23. Abramowitz, Milton, and Irene A. Stegun, Handbook of Mathematical Functions with Formulas,
Graphs, and Mathematical Tables, New York, NY, Dover Publications, Inc., 1965, Chapter 26.
24. Mortari, Daniele, “ESOQ2 Single-Point Algorithm for Fast Optimal Attitude Determination,”
Paper AAS 97–167, AAS/AIAA Space Flight Mechanics Meeting, Huntsville, AL, February 10-12, 1997.
25. Mortari, Daniele, “Second Estimator of the Optimal Quaternion,” Journal of Guidance, Control, and
Dynamics, Vol. 23, No. 5, September-October 2000, pp. 885-888.
26. The Math Works, Inc., MATLAB User’s Guide, Natick, MA, 1992.
27. Deutschmann, J., “Comparison of Single Frame Attitude Determination Methods,” Goddard Space
Flight Center Memo to Thomas H. Stengle, July 26, 1993.

View publication stats