FAST MATRIX BASED COMPUTATION OF EIGENVALUES IN POLSAR DATA

We describe calculation of eigenvalues of 2×2 and 3×3 Hermitian matrices as used in the analysis of multilook polarimetric SAR data. The eigenvalues are calculated as the roots of quadratic or cubic equations. The methods are well suited for fast matrix oriented computer implementation and the speed-up over calculations based on for-loops and built-in eigenproblem solvers is enormous.


INTRODUCTION
This paper deals with fast matrix based calculations of eigenvalues of 2×2 and 3×3 Hermitian matrices. Specifically, we compare (Matlab) computing times for the suggested matrix oriented calculations with a simple implementation with for loops over rows and columns and calls to the built-in eigenvalue function eig.
We use a change detection setting for illustration of the calculations. In the analysis of multi-look polarimetric synthetic aperture radar (SAR) data described by Hermitian (variance-) covariance matrices, the complex Wishart distribution can be used for change detection between acquisitions at two time points [1]. Many authors have worked with this type of change detection, see for example [2][3][4][5][6][7]. In [3,5] we work specifically with change between several time points.
Used in combination with a polarimetric change detector the Loewner order [8] calculates whether the difference between two Hermitian covariance matrices is positive definite, negative definite or indefinite. This can be done by means of eigenvalues or pivots [9,10]. The Loewner order allows us to establish whether the radar response increases, decreases or if it changes structure or nature between two time points in a way such that we cannot say whether it has increased or decreased. The Loewner order is thus a measure of direction of detected change.
The data used here to illustrate the techniques come from the Danish airborne EMISAR system [11,12].

MULTILOOK POLARIMETRIC SAR DATA
In the covariance matrix formulation of multilook polarimetric SAR image data each pixel can be described by a complex 3×3 matrix This matrix is Hermitian also known as self-adjoint, i.e., the matrix is equal to its own conjugate transpose. If we multiply by the number of looks, n, Z = n C will follow a complex Wishart distribution (for fully developed speckle). This is the matrix that is used in the change * alan@dtu.dk, https://people.compute.dtu.dk/alan detection methods described in [2,3,5,7]. Below and for the methods in [9,10] we may work on C or Z, here we'll use Z.

EIGENVALUES
Let us use the notation The first, fourth, fifth and sixth terms are real. The second and third terms are complex conjugates of each other. Since the imaginary parts cancel out, the sum of the two terms is twice the real part. Both trace and determinant of Z are real. Fast matrix inversion and fast determinant calculation are dealt with in [13].
Below we give computational details for determination of the eigenvalues of Z. For completeness we start by giving the more well known solution to the dual pol case. The solutions given are well suited for fast matrix based computer implementation.

Dual pol case
For dual polarimetry we have, say, the upper left 2×2 matrix of Z only k a a * ξ .
The trace is k + ξ, the determinant is kξ − |a| 2 , both are real. The eigenvalue problem can be written as which has the well known solution Here A = 1, B is minus the trace and C the determinant so we get The discriminant (k − ξ) 2 + 4|a| 2 is always nonnegative so the eigenvalues are real. The discriminant is zero if and only if k = ξ and |a| 2 = 0, i.e., the matrix we are analyzing is proportional to the identity matrix in which case the double eigenvalue is k = ξ.

Quad/full pol case
In this case we have the full matrix Z. The eigenvalues λ of Z are the roots of the cubic equation The eigenvalues are important for example in the Cloude-Pottier decomposition of polarimetric SAR data [14].
To find the inflection point (λi, f (λi)) of the cubic function with first and second derivatives set the second derivative to zero leading to To solve the cubic equation f (λ) = 0, divide by A and introduce which translates the inflection point along the x-axis so that its abscissa becomes zero to obtain is a so-called depressed cubic (which has no quadratic term) where At the inflection point, which translates the inflection point along the y-axis so that its ordinate becomes zero to see that g(−x) = −g(x). This shows that the cubic function is point symmetric (i.e., it has 180 • rotational symmetry) around the inflection point.
Inserting our A, B and C into the expression for 3p we get Since all terms are negative squares, p ≤ 0 and p = 0 if and only if k = ξ = ζ and |a| 2 = |b| 2 = |ρ| 2 = 0, i.e., the matrix we are analyzing is proportional to the identity matrix in which case the triple eigenvalue is k = ξ = ζ.
For the critical points xc of f (x) we have f (xc) = 3x 2 c + 3p = 0, so xc = ± √ −p. p < 0 means we have an inflection point and both a local minimum (at xc = √ −p) and a local maximum (at xc = − √ −p). For p ≥ 0 we have an inflection point only, no local extrema.
To find the eigenvalues substitute (due to François Viète, late 1500s) into x 3 + 3px + 2q = 0 to get which leads to a quadratic equation for z 3 The solution for z 3 is The solution depends on the sign of the discriminant q 2 + p 3 (q 2 + p 3 < 0 is referred to as the casus irreducibilis, the irreducible case). Instead we may use another substitution (also due to Viète) for x in x 3 + 3px + 2q = 0 namely x = u cos θ leading to u 3 cos 3 θ + 3pu cos θ + 2q = 0.
Under the heading "Cubic equation", Wikipedia has a good description of the problem and a great illustration of the above trigonometric solution to the cubic equation. This illustration is reproduced in Figure 1.
Note, that with the cosine substitution we easily see that for p < 0 the eigenvalues are real (non-complex) irrespective of the sign of the discriminant q 2 + p 3 .

Casus irreducibilis
In the casus irreducibilis or the irreducible case, the discriminant q 2 + p 3 < 0 and we may obtain a solution without the above cosine substitution, [15]. We have z 3 may be written in polar coordinates with magnitude ρ 3 and argument 3θ The cube roots of z 3 k are we get This is the same solution as obtained above with the cosine substitution. However, here we need p 2 + q 3 < 0 (which we have not checked), with the cosine substitution we need p < 0 only (which we have checked).

DATA EXAMPLE
The method is illustrated in a change detection setting as described in [1] with airborne EMISAR [11,12] data from 17 April and 20 May covering an agricultural test site near Foulum, Denmark, see Figure 2. The gray areas in the two images in the top row are wooded regions which, as opposed to most of the surrounding agricultural fields, show no change. Specifically, the bottom-right image shows significant change between the two time points colour coded for direction of change as determined by the Loewner order, [9,10].
Change and direction of change can be determined at patch or segment level also [5].

CONCLUSIONS
The speed-up factors (all approximate and for a few tests on 1024 by 1024 pixel images) for fast matrix based computer implementation based on the above eigenvalue calculations compared with a simple implementation based on calls to Matlab's built-in eigenvalue solver eig in for loops over rows and columns, is 350 for dual pol and 175 for quad/full pol, both enormous speed-ups.
The largest absolute value of the difference between the eigenvalues obtained from the two methods is less than 10 −11 .
Matlab code covering quad/full pol as well as azimuthal symmetry and dual pol including diagonal only data for change detection in polSAR data (with support functions), for calculating eigenvalues and derived quantities (entropy and two anisotropies) is available on the author's homepage. For 1024 by 1024 pixels full/quad pol, the eigenvalue calculations (with auxiliary variables entropy and two anisotropies) take a little less than 0.12 seconds elapsed time carried Polarimetric airborne EMISAR data (13 looks) for the two time points 17 April (X) and 20 May (Y ), 1024 rows by 1024 columns 5m pixels, somewhat untraditionally we show the logarithms of the descending eigenvalues as RGB (top row, where the gray areas in the two images are wooded regions which exhibit no change), the test statistic (−2ρ ln Q, high values, i.e., bright pixels indicate change, low values, i.e., dark pixels indicate no change) for change detected between the two time points as described in [1] (middle row left), the associated p-value (the change probability) thresholded at 99% (middle row right), positive definite matrix difference in red (i.e., Y <L X so X dominates, see [9,10]), negative definite matrix difference in green (i.e., X <L Y so Y dominates), indefinite matrix difference in yellow (bottom row left), and the same combination where the p-values in the Wishart based test are larger than 99% (bottom row right). The background gray scale image in the no-change regions is the temporal mean of the leading eigenvalues at the two time points. out with Matlab R2019b on a MacBook Pro from 2019, 2.3 GHz 8-Core Intel i9 processor, 64 GB 2667 MHz DDR4 memory.
The methods described here can be used in the analysis of Sentinel-1 and Radarsat-2 data also, as well as in other contexts, for example in the analysis of real symmetric variance-covariance matrices from RGB imagery.