1

I have an algorithm for calibrating a vector magnetometer. The input is $N$ readings of the $x$, $y$, $z$ axes: $(x_1, x_2, \dots, x_N)$, $(y_1, y_2, \dots, y_n)$, and $(z_1, z_2, \dots, z_N)$. The algorithm fits an ellipsoid to the data by estimating a symmetric $3\times3$ matrix $A$. In order to calibrate the system, it needs to calculate $\sqrt{A}$. I am adapting the algorithm for a microcontroller with very little memory, so cannot load standard matrix manipulation libraries. Is there an explicit formula for calculating the square root of $3\times3$ positive definite matrix?

Mohsen Shahriari
  • 6,718
Obromios
  • 121
  • 5

2 Answers2

3

There may not be an explicit formula, but there may be some methods that give a fair approximation to the square root.

For instance, if you know that the eigenvalues of $A$ are in $[0,1]$, you can take a polynomial approximation $P(x)$ to $\sqrt{x}$ on $[0,1]$, and then consider $P(A)$ that approximates $\sqrt{A}$. We can always reduce to "eigenvalues of $A$ in $[0,1]$" by multiplying by some positive constant.

Another way is to use the Newton approximation to the solution of the equation $x^2 = A$. Start with an initial term say $x_0 = I$, $x_{n+1} = \frac{1}{2}( x_n+ A \cdot x_n^{-1})$, it should converge quickly to $\sqrt{A}$.

orangeskid
  • 53,909
  • This is most useful thank you, and I may end up using it as my approach. I have upvoted your answer, but have accepted user1551's answer because it does match exactly the question I asked. – Obromios May 16 '21 at 21:13
  • @Obromios: You are very welcome! The second method should provide a very fast convergence. Basically, once the sequence stabilizes, you probably hit $\sqrt{A}$. It works for positive matrices $A$ in any dimensions, and also for other equations, say like $x^m = A$, $x$ positive. – orangeskid May 16 '21 at 21:38
  • I have added more information to the answer, which makes it is clearer that your eigenvalue approach could work as well. The maximum eigenvalue will the longest principle axis, which will be given by the reading with the maximum $sqrt{x^2_i + y^2_i + z^_i}$, so dividing by this maximum magnitude should ensure the eigenvallues are on [0,1]. – Obromios May 17 '21 at 00:23
  • @Obromios: Just found this post, > 20 years old... https://dspguru.com/dsp/tricks/square-root/ looks like the same approaches. That is a confirmation.... Now, about a polynomial approximation to the square root function. The Taylor series truncation for $\sqrt{1-x}$ is not good at $x=1$. So perhaps use other ways to get a poly approximation. Still I think that the recurrence would work very fast, get the result up to the error in calculations. – orangeskid May 17 '21 at 00:53
  • Just tried the Newton method, amazingly quick convergence. I think I will go with that, but checking for any problems on real data. Thank you – Obromios May 17 '21 at 01:01
  • If you think the question might be of use to others, you might like to upvote it. – Obromios May 17 '21 at 01:02
  • 1
    @Obromios: I did, but it seems somebody from the closing brigade gave it a -1. Glad that Newton works! – orangeskid May 17 '21 at 01:09
3

I would rather use an iterative method (e.g. QR iteration or Jacobi's algorithm) to diagonalise $A$ and find its square root, but there is a semi-closed-form formula for $B=\sqrt{A}$. By Cayley-Hamilton theorem, $$ B^3-\operatorname{tr}(B)B^2+\operatorname{tr}(B^{-1})\det(B)B-\det(B)I=0.\tag{1} $$ Multiply both sides of $(1)$ by $B$, we obtain $$ B^4-\operatorname{tr}(B)B^3+\operatorname{tr}(B^{-1})\det(B)B^2-\det(B)B=0.\tag{2} $$ Substitute $(1)$ into $(2)$ and using the fact that $B^2=A$, we obtain $$ A^2+\operatorname{tr}(B)\left[-\operatorname{tr}(B)A+\operatorname{tr}(B^{-1})\det(B)B-\det(B)I\right]+\operatorname{tr}(B^{-1})\det(B)A-\det(B)B=0. $$ Therefore \begin{align} B&=\frac{-A^2+\left(\operatorname{tr}(B)^2-\operatorname{tr}(B^{-1})\det(B)\right)A+\operatorname{tr}(B)\det(B)I} {\left(\operatorname{tr}(B)\operatorname{tr}(B^{-1})-1\right)\det(B)}\\ &=\frac{-A^2+(\alpha+\beta)A+\gamma\delta I} {\alpha\gamma-\delta}\tag{3}\\ \end{align} where \begin{align} \alpha&=\sqrt{\lambda_1\lambda_2}+\sqrt{\lambda_2\lambda_3}+\sqrt{\lambda_3\lambda_1},\\ \beta&=\operatorname{tr}(A),\\ \gamma&=\sqrt{\lambda_1}+\sqrt{\lambda_2}+\sqrt{\lambda_3},\\ \delta&=\sqrt{\det(A)} \end{align} and $\lambda_1,\lambda_2,\lambda_3$ are the eigenvalues of $A$.

So, to find $B$ using $(3)$, we need the eigenvalues $A$. Let $s=\operatorname{tr}(A^2)$. The characteristic equation of $A$ is then $$ x^3-\beta x^2+\frac12(\beta^2-s)x-\delta^2=0\tag{4} $$ and $\lambda_1,\lambda_2,\lambda_3$ are its roots. Now $(4)$ can be solved by radicals.

user1551
  • 139,064