Wikipedia:Reference desk/Archives/Mathematics/2012 August 7

Welcome to the Wikipedia Mathematics Reference Desk Archives
The page you are currently viewing is an archive page. While you can leave answers for any questions shown below, please ask new questions on one of the current reference desk pages.

August 7

Eigencomposition of real, asymmetric tridiagonal matrices

My linear algebra is rusty... Is there a simple and efficient way to find the eigendecomposition

$\mathbf {A} =\mathbf {Q} \cdot {\boldsymbol {\Lambda }}\cdot \mathbf {Q} ^{-1}$

of a n×n real-valued, asymmetric, diagonalizable, tridiagonal matrix A? All columns in A sum up to zero. All diagonal elements are negative, and all sub- and superdiagonal elements are positive. Thus, I know that the largest eigenvalue $\lambda _{1}\equiv 0$. n is in the range order 5-30. A typical example of A could be this matrix

[[-0.0733  0.0544  0.      0.      0.      0.      0.      0.      0.      0.    ]
 [ 0.0733 -0.1226  0.0068  0.      0.      0.      0.      0.      0.      0.    ]
 [ 0.      0.0682 -0.0931  0.0568  0.      0.      0.      0.      0.      0.    ]
 [ 0.      0.      0.0863 -0.1185  0.0251  0.      0.      0.      0.      0.    ]
 [ 0.      0.      0.      0.0617 -0.1658  0.1971  0.      0.      0.      0.    ]
 [ 0.      0.      0.      0.      0.1407 -0.264   0.0093  0.      0.      0.    ]
 [ 0.      0.      0.      0.      0.      0.0669 -0.03    0.0069  0.      0.    ]
 [ 0.      0.      0.      0.      0.      0.      0.0207 -0.0117  0.0029  0.    ]
 [ 0.      0.      0.      0.      0.      0.      0.      0.0048 -0.0095  0.0511]
 [ 0.      0.      0.      0.      0.      0.      0.      0.      0.0066 -0.0511]]

I have found a simple LU decomposition method for tridiagonal matrices

$\mathbf {A} =\mathbf {L} \cdot \mathbf {U}$

by which A can be decomposed into two tridiagonal matrices, where the superdiagonal elements are zero in L and the subdiagonal elements are zero in U and all diagonal elements are unity. Could this be of use?

I also found a §11.3 "Eigenvalues and Eigenvectors of a Tridiagonal Matrix" in Numerical recipes. However, §11.3 only seems to cover symmetric matrices.

I also have a copy of Matrix computations by Golub and van Loan, but I did not find any coverage of the topic there either. --Slaunger (talk) 14:03, 7 August 2012 (UTC)

Rotational transformation on a map

I'm working on a project with Google Maps, which uses a Mercator projection.

Imagine I have one point on the map at position lat1, lng1.

I have another point on the map at position lat2, lng2.

What I want to do is find lat3, lng3, which is defined as being the same distance from lat1, lng1 as lat2, lng2 is, but rotate an arbitrary angle (say, -16 degrees).

It should look something like this, but on a map:

|
|           3
|
|              2
| 
|
|      1
|______________________

How should I go about doing this? Obviously this is a rotational transformation on a 3-D surface (the Earth, and the distances for this application are sufficient than a 2-D abstraction will not work correctly).

If you could explain this in code or pseudocode it would be most helpful, as opposed to plunking down a matrix transformation or linking to this article which I have read but don't quite understand how to apply. The actual code is written in Javascript, which does not have any kind of standard rotational matrix function, so please don't rely on one of those...

Thanks in advance. --Mr.98 (talk) 20:32, 7 August 2012 (UTC)

See spherical trigonometry. Bo Jacoby (talk) 21:16, 7 August 2012 (UTC).

Thanks. But can you give me some more practical pointers? The article is quite theoretical. --Mr.98 (talk) 21:33, 7 August 2012 (UTC)

You don't need spherical geometry. Simply convert them all to unit vectors then use Rodrigues' rotation formula. If you wanted you could generate the rotation (as a matrix or quaternion) for the rotation then use that, but for a single rotation this formula should work.--JohnBlackburne^words_deeds 01:08, 8 August 2012 (UTC)

Cool. Can you help me turn this into a function? Pseudocode? Anything other than just a page full of equations. --Mr.98 (talk) 01:27, 8 August 2012 (UTC)

Can you obtain the vector, at point 1, normal to the sphere ? This will be your axis of rotation. If you know the center point of the sphere, this should be simple, as the vector is simply the line connecting the center point and point 1. StuRat (talk) 01:34, 8 August 2012 (UTC)

Call the points A,B,C rather than 1,2,3. Let P be the north pole.

First consider the triangle PAB having angles P₁=lng2−lng1, A₁, B₁ opposite the sides p₁, a₁=90°−lat2, b₁=90°−lat1.

The spherical law of cosines gives p₁, and the law of sines#spherical case gives A₁

cos p₁ = cos a₁ cos b₁ + sin a₁ sin b₁ cos P₁

sin A₁/sin a₁ = sin P₁/sin p₁

Next consider triangle PAC having angles P₂, A₂=A₁−16°, C₂ opposite the sides p₂=p₁, a₂, c₂=b₁.

The law of cosines gives a₂, and the law of sines gives P₂

cos a₂ = cos p₂ cos c₂ + sin p₂ sin c₂ cos A₂

sin P₂/sin p₂ = sin A₂/sin a₂

Then lat3=90°−a₂ and lng3=lng1+P₂ is your final answer. Bo Jacoby (talk) 07:20, 8 August 2012 (UTC).

Thanks; I'll see if I can make this work. And just to emphasize it again: I'm not a math guy. Wikipedia articles on math in and of themselves don't mean much to me, because they quickly drop into jargon I don't understand and formalisms I've never used before. I can make things work in code, because that's just plugging functions together in a way that I recognize, but beyond that, I'm pretty clueless. (And I should emphasize, perhaps, that these are not homework problems in any way, shape, or form!) So I do appreciate your efforts to make it slightly more programmatic for me! --Mr.98 (talk) 16:23, 8 August 2012 (UTC)

I'm sure Bo Jacoby's solution is mathematically correct but I'm worried about its stability—it may sometimes give wildly wrong answers in practice because of roundoff error. Here's my attempt at a solution that will never go too wrong, following JohnBlackburne's suggestion:

$\mathbf {k} =(\cos \mathrm {lat} _{1}\cos \mathrm {lng} _{1},\;\cos \mathrm {lat} _{1}\sin \mathrm {lng} _{1},\;\sin \mathrm {lat} _{1})$

$\mathbf {v} =(\cos \mathrm {lat} _{2}\cos \mathrm {lng} _{2},\;\cos \mathrm {lat} _{2}\sin \mathrm {lng} _{2},\;\sin \mathrm {lat} _{2})$

$\mathbf {w} =\mathbf {v} \cos \theta +(\mathbf {k} \times \mathbf {v} )\sin \theta +\mathbf {k} (\mathbf {k} \cdot \mathbf {v} )(1-\cos \theta )$

$\mathrm {lat} _{3}=\mathrm {atan2} (w_{z},{\sqrt {w_{x}^{2}+w_{y}^{2}}})$

$\mathrm {lng} _{3}=\mathrm {atan2} (w_{y},w_{x})$

Here's the computation of w in components:

$\mathbf {k} \cdot \mathbf {v} =k_{x}v_{x}+k_{y}v_{y}+k_{z}v_{z}$

$w_{x}=v_{x}\cos \theta +(k_{y}v_{z}-k_{z}v_{y})\sin \theta +k_{x}(\mathbf {k} \cdot \mathbf {v} )(1-\cos \theta )$

$w_{y}=v_{y}\cos \theta +(k_{z}v_{x}-k_{x}v_{z})\sin \theta +k_{y}(\mathbf {k} \cdot \mathbf {v} )(1-\cos \theta )$

$w_{z}=v_{z}\cos \theta +(k_{x}v_{y}-k_{y}v_{x})\sin \theta +k_{z}(\mathbf {k} \cdot \mathbf {v} )(1-\cos \theta )$

-- BenRG (talk) 18:41, 8 August 2012 (UTC)

With a numeric example (lat1, lng1, lat2, lng2, rotation angle) we can try both methods and compare results and clarify the suspicion of instability. Spherical trigonometry has proved useful to astronomers during more than 500 years. Bo Jacoby (talk) 20:37, 8 August 2012 (UTC).

Considering that the earth is better approximated by an ellipsoid than a sphere isn't the usability of spherical trigonometry tied up to if a spherical approximation is reasonable? The relevance of a spherical assumption depends on the precision required in the solution and the distances this should work for relative to the size of the earth. I know from own work with radar systems, where you need to convert a slant range and bearing measurement relative to a known lat and long for a sensor origin into a lat, long for the measurement, that this is nontrivial and involves also considering the eccentricity of the reference ellipsoid if it has to be really precise (within 1 meter of the correct position) and the slant range is a few houndred kilometers. --Slaunger (talk) 05:34, 10 August 2012 (UTC)

The solution with spherical trigonometry has problems when two of the points P, A, B are identical or antipodal. E.g. if A and B are antipodal then sin P₁/sin p₁ is 0/0. I think my solution has no such special cases. I'm a little worried about precision, though. An IEEE 754 double has 53 bits of mantissa, and R_e / 2⁵³ ≈ 1 nm, which is fine, but the rule of thumb is that you need half the bits for "support", and R_e / 2^53/2 ≈ 10 cm, which isn't so great. The full 64-bit mantissa of the x87 would be plenty, but the trend these days is to disable it in the name of reproducibility on other, inferior processors. -- BenRG (talk) 07:45, 10 August 2012 (UTC)

To Slaunger: The OP gave the example '-16 degrees' indicating than 0.1 degrees accuracy is enough. To BenRG: If A and B are antipodal and C has the same distance from A as B, then the solution is trivial: C=B. Bo Jacoby (talk) 11:57, 10 August 2012 (UTC).

Trivial cases such as this generally need special handling (and hence code) in software, and should not be so glibly discounted. Handling of instability (including identifying it) near such cases frequently needs expertise, and will often be neglected due to lack of awareness thereof on the part of the programmer. Formulations that are stable throughout the domain of the problem are thus to be preferred where these are available. Consider the embarassment of the NASA upon ignoring the problem of gimbal lock; Google's use might be less critical but unnecessary erratic behaviour may retain an irritating presence for ages. A formulation that maps the problem into 3D space with Cartesian coordinates provides global stability in a straightforward manner with no special cases. BenRG's suggestion is such a formulation (although the final translation to calculate longitude, not stated by BenRG, is unavoidably unstable near the north and south poles in a relatively benign way but should nevertheless be handled to avoid failure). — Quondum^☏ 17:28, 10 August 2012 (UTC)

Perfect is the enemy of good. Bo Jacoby (talk) 21:17, 10 August 2012 (UTC).

I don't understand what you mean by linking that article, but what Quondum said is true: calculations that have singularities over the (exact mathematical) reals behave badly in some neighborhood of the singularity when using fixed-precision floats. It's possible to use a patchwork of two or more algorithms but it's very hard to pick good boundaries and avoid discontinuities caused by different rounding errors at the boundaries.

Quondum, I did overlook that, but Mr.98 is using Javascript, and at least in standard ECMAscript Math.atan2 is defined to always return an angle unless an argument is NaN. In C atan2(0,0) is allowed to trap or return NaN, so you would need to check for that case and return a dummy angle instead. All other argument values are fine because atan2 was written by experts (one hopes). -- BenRG (talk) 01:47, 11 August 2012 (UTC)

I mean that the solution using spherical trigonometry is good. The OP may use it and find correct results. Theoretically singularities occur, but with zero probability. Spending attention to such situations in order to make the method perfect is derailing the discussion. So Perfect is the enemy of Good. Bo Jacoby (talk) 20:52, 11 August 2012 (UTC).

Hmm, instead of the above maybe you should have a look at this online calculator, source code, and algorithms in linked pdf.—eric 18:30, 10 August 2012 (UTC)