OSA's Digital Library

Optics Express

Optics Express

  • Editor: C. Martijn de Sterke
  • Vol. 20, Iss. 1 — Jan. 2, 2012
  • pp: 17–31
« Show journal navigation

Mueller matrix roots algorithm and computational considerations

H. D. Noble and R. A. Chipman  »View Author Affiliations


Optics Express, Vol. 20, Issue 1, pp. 17-31 (2012)
http://dx.doi.org/10.1364/OE.20.000017


View Full Text Article

Acrobat PDF (1331 KB)





Browse Journals / Lookup Meetings

Browse by Journal and Year


   


Lookup Conference Papers

Close Browse Journals / Lookup Meetings

Article Tools

Share
Citations

Abstract

Recently, an order-independent Mueller matrix decomposition was proposed in an effort to elucidate the nine depolarization degrees of freedom [Handbook of Optics, Vol. 1 of Mueller Matrices (2009)]. This paper addresses the critical computational issues involved in applying this Mueller matrix roots decomposition, along with a review of the principal matrix root and common methods for its calculation. The calculation of the pth matrix root is optimized around p = 105 for a 53 digit binary double precision calculation. A matrix roots algorithm is provided which incorporates these computational results. It is applied to a statistically significant number of randomly generated physical Mueller matrices in order to gain insight on the typical ranges of the depolarizing Matrix roots parameters. Computational techniques are proposed which allow singular Mueller matrices and Mueller matrices with a half-wave of retardance to be evaluated with the matrix roots decomposition.

© 2011 OSA

1. Introduction

Polarization elements and their associated Mueller matrices are typically decomposed and analyzed in terms of the polarization properties of diattenuation, retardance, and depolarization. For Mueller matrices with a mixture of all three, the properties are distributed among the matrix elements in a complex manner. Lu and Chipman described a Mueller matrix data reduction method based on polar decomposition [1

1. S.-Y. Lu and R. A. Chipman, “Interpretation of mueller matrices based on polar decomposition,” J. Opt. Soc. Am. A 13, 1106–1113 (1996). [CrossRef]

]. It decomposes the Mueller matrix into an order-dependent product of a depolarizer, a retarder, and a diattenuator. Ossikovski’s symmetric decomposition method decomposes a depolarizing Mueller matrix into a sequence of five factors: a diagonal depolarizer between two retarder and diattenuator pairs [2

2. R. Ossikovski, “Analysis of depolarizing mueller matrices through a symmetric decomposition,” J. Opt. Soc. Am. A 26, 1109–1118 (2009). [CrossRef]

]. The diattenuation, retardance, and depolarization parameters calculated by these methods depend on the arbitrary order of the decomposition. The Cloude additive decomposition [3

3. S. R. Cloude, “Conditions for the physical realisability of matrix operators in polarimetry,” Proc. SPIE 1166, 177–185 (1989).

] separates a Mueller matrix into the sum of four non-depolarizing Mueller matrices which are scaled by eigenvalues of the Mueller matrix covariance matrix. Since this is an additive decomposition, its terms are order-independent. None of these methods clearly elucidate the nine degrees of freedom associated with depolarization.

Recently, a Mueller matrix roots decomposition was introduced by Chipman in [4

4. R. A. Chipman, Handbook of Optics, Vol. 1 of Mueller Matrices (McGraw Hill, 2009), 3rd ed.

] to provide an order-independent description of polarization properties and provide a clear analysis of the nine depolarization degrees of freedom. The Mueller matrix roots decomposition extends the concept of the Jones N-matrices [5

5. R. C. Jones, “A new calculus for the treatment of optical systems. vii. properties of the n-matrices,” J. Opt. Soc. Am. 38, 671–683 (1948). [CrossRef]

] to Mueller matrices by dividing the Mueller matrix M into p infinitesimal slices:
N=Mp,
(1)
as illustrated in Fig. 1, where N is the Mueller matrix for one small slice. When p becomes very large, the polarization properties of Mp separate (and become order-independent) as the diattenuation, retardance, and depolarization become very small near the identity matrix. In the limit as p approaches infinity, the principal pth root of a large class of Mueller matrices approaches the identity matrix I,
limpMp=I.
(2)
In this study, the class of Mueller matrices that obey Eq. (2) are called uniform Mueller matrices. The Mueller matrix roots decomposition from [4

4. R. A. Chipman, Handbook of Optics, Vol. 1 of Mueller Matrices (McGraw Hill, 2009), 3rd ed.

] is applicable to this subset of Mueller matrices, with the exception of the non-uniform special cases highlighted in section 4.

Fig. 1 Taking the root of a uniform Mueller matrix is analogous to slicing it into very thin identical pieces.

The calculation of the principal pth root of square matrices is an extensively studied subject [6

6. N. J. Higham, Functions of Matrices: Theory and Computation (Society for Industrial and Applied Mathematics, 2008). [CrossRef]

8

8. C.-H. Guo, “On newton’s method and halley’s method for the principal pth root of a matrix,” Linear Algebra Appl. 432, 1905–1922 (2010). [CrossRef]

], and numerical accuracy and noise are well understood within the field of numerical computing [9

9. R. D. Skeel and J. B. Keiper, Elementary Numerical Computing with Mathematica, Chapter 2 (McGraw Hill, 1993).

]. In this study, these concepts are applied to the calculation of the Mueller matrix roots for the purpose of Mueller matrix decomposition. In order for the Mueller matrix roots decomposition from [4

4. R. A. Chipman, Handbook of Optics, Vol. 1 of Mueller Matrices (McGraw Hill, 2009), 3rd ed.

] to be order-independent, the pth root must be optimized with all of these topics in mind. This paper reviews and updates the Mueller matrix roots decomposition from [4

4. R. A. Chipman, Handbook of Optics, Vol. 1 of Mueller Matrices (McGraw Hill, 2009), 3rd ed.

], addresses several computational issues, and highlights several common non-uniform special cases that arise. Considerations for these computational issues and special cases are implemented in a matrix roots decomposition algorithm, which is applied in a statistical analysis of a large quantity of randomly generated physical Mueller matrices.

2. Matrix roots decomposition

The goal of the Mueller matrix roots decomposition is to calculate the magnitudes of sixteen distinct properties of the Mueller matrix, including the nine depolarization degrees of freedom, for the set of uniform Mueller matrices.

To calculate the matrix roots decomposition of a Mueller matrix M, a Mueller matrix N with infinitesimal polarization properties is first calculated from the pth principal root of M, where p is some large integer (typically 105),
N=Mp.
(3)
The appropriate choice of p is discussed in detail in section 5. If p is too small, the Mueller matrix roots decomposition fails to be order-independent.

The infinitesimal polarization parameters d0 through d15 are defined from the symmetric and antisymmetric parts of N:
N=d0(1d1+d7d2+d8d3+d9d1d71f13d6+d12d5+d11d2d8d6+d121f14d4+d10d3d9d5+d11d4+d101f15),
(4)
where d13, d14, and d15 are solved from the first-order generator products in terms of the parameters f13, f14, and f15:
d13=f14f132
(5)
d14=23(2f15+f14+f13)
(6)
d15=16(f15+f14+f13).
(7)
The infinitesimal polarization parameters d0 through d15 are rescaled by p to produce the matrix roots parameters D0 through D15:
Dn=pdn,n=0,1,2,,15.
(8)
D0 through D15 parameterize the sixteen degrees of freedom of M.

There are three matrix roots parameters for diattenuation, D1, D2, D3, and three matrix roots parameters for retardance, D4, D5, D6. The three degrees of freedom for each property correspond to the axes in the Stokes/Mueller formalism (horizontal/vertical, 45°/135°, right/left circular).

Nine more parameters, D7 through D15, describe depolarizing effects - one for each of the nine depolarizing degrees of freedom described in [4

4. R. A. Chipman, Handbook of Optics, Vol. 1 of Mueller Matrices (McGraw Hill, 2009), 3rd ed.

]. There are three families of depolarizing parameters, each similarly divided into horizontal/vertical, 45°/135°, and right/left circular components. The terms for amplitude depolarization, D7, D8, and D9, share the same matrix elements as the parts of the Mueller matrix associated with diattenuation, on the horizontal/vertical, 45°/135°, and right/left circular axes [4

4. R. A. Chipman, Handbook of Optics, Vol. 1 of Mueller Matrices (McGraw Hill, 2009), 3rd ed.

]. They are named amplitude depolarization because they depolarize and affect the flux of an incident Stokes vector. The terms for phase depolarization, D10, D11, and D12, correspond with the parts of the Mueller matrix associated with retardance on the horizontal/vertical, 45°/135°, and right/left circular axes. They do not affect the flux of an incident Stokes vector. D13, D14, and D15 are named diagonal depolarization because they lie on the matrix diagonal. D15 expresses the overall isotropic depolarizing power, since D15 reduces the degree of polarization of all incident Stokes vectors equally, independent of the Stokes vector’s location on the Poincare sphere. D13 expresses the relative strength between the diagonal depolarization on the two linear axes (horizontal/vertical and 45°/135°). D14 expresses the relative strength between linear and circular diagonal depolarization. The diagonal depolarization parameters have been modified from [4

4. R. A. Chipman, Handbook of Optics, Vol. 1 of Mueller Matrices (McGraw Hill, 2009), 3rd ed.

] so that only the diagonal depolarization parameter D15 changes the Gil-Bernabeu depolarization index [10

10. R. A. Chipman, “Depolarization index and the average degree of polarization,” Appl. Opt. 44, 2490–2495 (2005). [CrossRef] [PubMed]

,11

11. J. J. Gil and E. Bernabeu, “Depolarization and polarization indices of an optical system,” Opt. Acta 33, 185–189 (1986). [CrossRef]

], while D0 through D14 do not [12

12. H. D. Noble, “Mueller matrix roots,” Doctoral dissertation (2011).

]. Table 1 lists the parameters D1 through D15 and categorizes them into their corresponding families and axes.

Table 1. The sixteen polarization properties of the Mueller matrix given by the Mueller matrix roots decomposition.

table-icon
View This Table

3. pth root

This section discusses the definition of the principal pth matrix root along with relevant examples and common methods of calculating the pth matrix root. This topic is of critical importance for calculating the Matrix roots, since matrices have multiple roots.

Matrices have multiple roots. For a nonsingular matrix A ∈ ℂn×n (in complex space) with s distinct eigenvalues, there are precisely ps pth roots [6

6. N. J. Higham, Functions of Matrices: Theory and Computation (Society for Industrial and Applied Mathematics, 2008). [CrossRef]

]. So long as the matrix A ∈ ℂn×n has no negative, real eigenvalues, there is a unique pth root of A whose eigenvalues’ arguments lie between −π/p and π/p, and that unique root is defined as the principal root of A [6

6. N. J. Higham, Functions of Matrices: Theory and Computation (Society for Industrial and Applied Mathematics, 2008). [CrossRef]

]. If A is real, then its principal root A1/p is real.

Singular matrices (such as polarizer Mueller matrices) do not have principal matrix roots. Nonetheless, matrix roots of singular Mueller matrices near the identity matrix can still be found. Methods for calculating roots of singular Mueller matrices (such as linear polarizers) are discussed in section 4.1. Methods for calculating roots of non-depolarizing and depolarizing Mueller matrices with π retardance are discussed in section 4.2.

3.1. Principal matrix root algorithms

This section provides a brief summary of common methods for calculating principal matrix roots, including the Schur method, Newton’s method, and the Schur-Newton algorithm. This provides a starting point for the reader who is interested in applying these methods to the calculation of Mueller matrix roots.

Schur methods form a Schur decomposition of A and compute a pth root of the resulting upper triangular factor using various (stable) recursive formulae [13

13. N. J. Higham, “Stable iterations for the matrix square root,” Num. Algor. 15, 227–242 (1997). [CrossRef]

]. Newton’s method calculates the pth root of A using an iterative approach [6

6. N. J. Higham, Functions of Matrices: Theory and Computation (Society for Industrial and Applied Mathematics, 2008). [CrossRef]

]. The Newton method is largely considered to have poor convergence and stability properties [14

14. M. H. Smith, “Optimization of a dual-rotating-retarder mueller matrix imaging polarimeter,” Appl. Opt. 41, 2488–2493 (2002). [CrossRef] [PubMed]

], and the Schur method from [14

14. M. H. Smith, “Optimization of a dual-rotating-retarder mueller matrix imaging polarimeter,” Appl. Opt. 41, 2488–2493 (2002). [CrossRef] [PubMed]

] is the numerically stable benchmark against which other methods are often compared [7

7. D. A. Bini, N. J. Higham, and B. Meini, “Algorithms for the matrix pth root,” Num. Algor. 39, 349–378 (2005). [CrossRef]

]. The Schur-Newton algorithm applies iterative computations to the upper triangular matrices from the Schur decomposition [6

6. N. J. Higham, Functions of Matrices: Theory and Computation (Society for Industrial and Applied Mathematics, 2008). [CrossRef]

]. Many of these algorithms are available in the MATLAB Matrix Computation Toolbox [7

7. D. A. Bini, N. J. Higham, and B. Meini, “Algorithms for the matrix pth root,” Num. Algor. 39, 349–378 (2005). [CrossRef]

].

Mathematica has built-in routines that diagonalize a matrix to easily calculate its root, so long as the matrix is diagonalizable. The diagonalization factors the matrix A into the product of three matrices,
A=Zdiag(λi)Z1,
(9)
where λi are the eigenvalues of A, the columns of Z are its eigenvectors, and diag(λi) is a diagonal matrix with its ith diagonal element equal to λi. Then the matrix root of A is calculated as
A1/p=Zdiag(λi)1/pZ1.
(10)
This method often (but not always) yields the principal root of a diagonalizable Mueller matrix, so long as its principal root exists. However if A is real and has some complex eigenvalues, then the computed A1/p, which should be real, may acquire a tiny imaginary part due to computational rounding errors. This imaginary part should be discarded. However, numerical instability can produce a large spurious imaginary part, so diagonalization-based computations of matrices with any imaginary eigenvalues should be treated with care [6

6. N. J. Higham, Functions of Matrices: Theory and Computation (Society for Industrial and Applied Mathematics, 2008). [CrossRef]

].

4. Special cases

4.1. Polarizers

The Mueller matrix for a polarizer is singular, regardless of ellipticity or orientation. Therefore it is not uniform, and the procedure described in section 2 will not yield meaningful matrix roots parameters when applied to a polarizer Mueller matrix. However, the treatment proposed in this section allows for this special case to be analyzed using the Mueller matrix roots decomposition by replacing the polarizer with an imperfect polarizer.

Because an ideal homogeneous polarizer [15

15. S.-Y. Lu and R. A. Chipman, “Homogeneous and inhomogeneous jones matrices,” J. Opt. Soc. Am. A 11, 766–773 (1994). [CrossRef]

] is always idempotent, the matrix roots decomposition will break down at the point where the pth root is calculated. For example, the formula for a homogeneous linear polarizer as a function of orientation θ is
LP(θ)=12(1cos2θsin2θ0cos2θcos22θ12sin4θ0sin2θ12sin4θsin22θ00000).
(11)
Note that
LP(θ)LP(θ)=LP(θ),
(12)
and thus
LP(θ)1/p=LP(θ).
(13)
The linear polarizer Mueller matrix is its own root to all orders, and the roots do not approach the identity matrix.

A homogeneous polarizer has two orthogonal, physical eigenpolarizations with eigenvalues (Tmax, 0), where Tmax is the maximum transmission, and the minimum transmission Tmin is 0. The polarizer can be replaced with a nearby uniform Mueller matrix by adjusting the maximum and minimum transmission by a small number ɛ. The partial polarizer matrix is now uniform, since the partial polarizer is a diattenuator with the same orthogonal eigenpolarizations, but the eigenvalues associated with its physical eigenpolarizations are (Tmaxɛ,ɛ).

This is accomplished by calculating the latitude η and orientation θ of the state of maximum transmission on the Poincare sphere from the polarizer’s three dimensional diattenuation vector {dH, d45, dR} [1

1. S.-Y. Lu and R. A. Chipman, “Interpretation of mueller matrices based on polar decomposition,” J. Opt. Soc. Am. A 13, 1106–1113 (1996). [CrossRef]

, 4

4. R. A. Chipman, Handbook of Optics, Vol. 1 of Mueller Matrices (McGraw Hill, 2009), 3rd ed.

], where
θ=12arctandHd45,and
(14)
η=arcsindRdH2+d452+dR2.
(15)
η and θ are then plugged in to the general formula for an elliptical diattenuator [4

4. R. A. Chipman, Handbook of Optics, Vol. 1 of Mueller Matrices (McGraw Hill, 2009), 3rd ed.

, 16

16. D. Goldstein, Polarized Light, 2nd ed. (Marcel Dekker, 2003). [CrossRef]

]. For an ideal polarizer, Tmax = 1 and Tmin = 0. The pth root of this partial polarizer approaches the identity matrix for large p - after replacing the polarizer with the nearby partial polarizer, the Mueller matrix roots decomposition algorithm can proceed as in the general case.

4.2. Half wave retarders

Mueller matrix roots for retarders are easily calculated and understood. For two linear retarders with the same fast-axis orientation, retardance is additive,
δ=δ1+δ2.
(16)
Thus the square root of a linear retarder is a retarder with half the retardance and the same fast-axis orientation. This is equivalent to cutting a wave plate in half. Similarly, the pth root of an ideal homogeneous linear retarder LR(δ,θ) with retardance δ at orientation θ is
LR(δ,θ)1/p=LR(δ/p,θ).
(17)
A retarder with a half-wave of retardance has negative, real eigenvalues, and therefore no principal root, so the half-wave retarder must be treated as a special case.

For example, the half-wave linear retarder oriented at 0°,
LR(π,0°)=(1000010000100001),
(18)
has negative real eigenvalues (λ = {−1,−1,1,1}), so the half wave retarder has no principal pth root. For any value of p, two eigenvalues of HWR1/p are (−1)p with arguments of −π/p and π/p, which lie on the edge of the principal root segment defined in section 3. While the desired solution in this case is not a principal root, it is the solution which approaches the identity matrix. The half wave retarder’s uniform pth root can be calculated analogously to the polarizer case, by means of a small perturbation.

Orientation θ and latitude η can be calculated from the three-dimensional retardance vector {δH, δ45, δR} [1

1. S.-Y. Lu and R. A. Chipman, “Interpretation of mueller matrices based on polar decomposition,” J. Opt. Soc. Am. A 13, 1106–1113 (1996). [CrossRef]

, 4

4. R. A. Chipman, Handbook of Optics, Vol. 1 of Mueller Matrices (McGraw Hill, 2009), 3rd ed.

] as follows,
θ=12arctanδHδ45,and
(19)
η=arcsinδRδH2+δ452+δR2.
(20)
The general form for an elliptical retarder with retardance δ, orientation θ, and latitude η can be written in terms of three linear retarders as
ER(η,θ,δ)=LR(η,θ+π/4)LR(δ,θ)LR(η,θπ/4),
(21)
where the general form for the linear retarder LR(δ,θ) with retardance δ and an orientation of θ is [4

4. R. A. Chipman, Handbook of Optics, Vol. 1 of Mueller Matrices (McGraw Hill, 2009), 3rd ed.

]
LR(δ,θ)=(10000cos22θ+cosδsin22θsin2δ2sin4θsinδsin2θ0sin2δ2sin4θcosδcos22θ+sin22θcos2θsinδ0sinδsin2θcos2θsinδcosδ).
(22)
After θ and η are calculated from the retardance vector, the half-wave retarder is perturbed by some small number (ɛ ≈ 10−7) to a nearby elliptical retarder ER′ of retardance (πɛ) by using equation 21:
ER(η,θ,δ=πɛ)=LR(η,θ+π/4)LR(πɛ,θ)LR(η,θπ/4).
(23)
This is a unitary transformation - therefore any higher order root of a half wave retarder can be calculated as follows:
ER(η,θ,δ=π)1p=LR(η,θ+π/4)LR((πɛ),θ)1pLR(η,θπ/4).
(24)
After following this procedure, the matrix roots of the perturbed half-wave retarders can be calculated without issue.

It is also possible to calculate the root in Eq. (24) without perturbation of the retardance,
HWR(η,θ)1p=LR(ɛ,θ+π/4)LR(πp,θ)LR(η,θπ/4),
(25)
but this method cannot be extended to a depolarizing Mueller matrix with a half-wave of retardance. This is discussed further in the following section.

4.3. Depolarizing non-uniform Mueller matrices

When half-wave retarders or polarizers are combined with other polarization properties, the resulting Mueller matrix is also non-uniform, and therefore cannot be analyzed with the standard Mueller matrix roots decomposition. Using the procedure discussed in this section, they can be perturbed to nearby uniform Mueller matrices.

A depolarizing Mueller matrix that also has a half-wave of retardance has negative, real eigenvalues and is therefore non-uniform. For example, the product of a half wave retarder oriented at 0° and a partial diagonal depolarizer with positive, real a, b, and c
(1000010000100001)(10000a0000b0000c)=(10000a0000b0000c),
(26)
has eigenvalues λ = {1,a,−b, −c}, and therefore no principal root. For cases where b = c, the perturbation approach used in section 4.2 can be modified to yield a uniform Mueller matrix. When bc, the perturbation approach fails to yield non-negative eigenvalues, and therefore the matrix will remain non-uniform.

To perturb a half-wave retarder mixed with depolarization and/or diattenuation properties (with equal negative eigenvalues) to a nearby uniform Mueller matrix, the Lu-Chipman decomposition [1

1. S.-Y. Lu and R. A. Chipman, “Interpretation of mueller matrices based on polar decomposition,” J. Opt. Soc. Am. A 13, 1106–1113 (1996). [CrossRef]

] is first calculated as an intermediate step,
M=MΔMRMD.
(27)
MΔ is a depolarizing Mueller matrix, MR is a pure retarder, and MD is a pure diattenuator. The half-wave retarder found in the MR term is perturbed to the nearby retarder MR with (πɛ) retardance as shown in equation 24. Then the terms are recombined to form M′, where
M=MΔMRMD.
(28)
The eigenvalues of M′ are no longer negative, since MΔ and MD always have positive eigenvalues [1

1. S.-Y. Lu and R. A. Chipman, “Interpretation of mueller matrices based on polar decomposition,” J. Opt. Soc. Am. A 13, 1106–1113 (1996). [CrossRef]

], and the eigenvalues of MR are positive. Therefore the uniform pth root of the perturbed depolarizing half-wave retarder from Eq. (28) can be calculated from M′, and its Mueller matrix roots decomposition parameters can be found.

Because the product of a polarizer and any other Mueller matrix is also singular, a polarizer multiplied by a depolarizer or retarder (or even a diattenuator) in any order is also singular. These Mueller matrices can be identified by testing for a diattenuation vector D of magnitude D = 1 and a depolarization index less than 1.

To replace a polarizer with a mixture of other polarization properties to a nearby uniform Mueller matrix, the Lu-Chipman decomposition [1

1. S.-Y. Lu and R. A. Chipman, “Interpretation of mueller matrices based on polar decomposition,” J. Opt. Soc. Am. A 13, 1106–1113 (1996). [CrossRef]

] is first calculated as an intermediate step,
M=MΔMRMD.
(29)
The polarizer’s maximum transmission can be adjusted by ɛ to the nearby uniform diattenuator (or partial polarizer) according to the procedure in section 4.1. The modified MD is substituted into Eq. (29) in place of MD, and the three matrices are recombined to form a modified, non-singular matrix M′:
M=MΔMRMD.
(30)
Then the pth root of M′ can be calculated, and its Mueller matrix roots decomposition parameters can be found.

5. Numerical accuracy and root order

5.1. Choice of p

The choice of root order p is an important consideration when calculating the matrix roots parameters. p must be large enough so that the Mueller matrix elements become sufficiently small for the polarization properties to separate and achieve and order-independent decomposition. However, the magnitude of p should also be as small as possible so as to minimize unnecessary loss of numerical precision from the calculation.

5.1.1. Accuracy

Because the Mueller matrix roots of retarders are straightforward and well-understood, they provide a good reference from which to study the convergence of Eq. (8). The pth root of a retarder Mueller matrix with retardance δ results in a pth principal matrix root with retardance δ/p. Therefore a Mueller matrix of an ideal elliptical retarder with retarder vector {δH,δ45,δR} should have matrix roots parameters of D4 = δH, D5 = δ45, and D6 = δR. In order to evaluate the correct choice of p in consideration of this criteria, the matrix roots retardance parameters (D4 through D6) were calculated (using Mathematica’s default double-precision machine arithmetic) for a pure elliptical retarder with retardance parameters {δH,δ45,δR} = {0.294,0.302,0.997} for different values of root order p. The relative error Δx was then calculated according to the following expression,
Δx=(δHD4)2+(δ45D5)2+(δRD6)2δ.
(31)

Figure 2 shows the relative error between the input retardance vector {δH,δ45,δR} and the corresponding matrix roots retardance vector D4 through D6 for different choices of p for the pth root. In Fig. 2, the relative retardance error converges to a minimum value just beyond the 105th root. For large p, the relative error increases, due to numerical rounding and the loss of precision associated with machine arithmetic.

Fig. 2 The error of the root calculation vs. root order for a simple retarder converges to a minimum relative error just beyond the 105th root.

An equally important factor to consider is the convergence of the Mueller matrix root parameter values. When the choice of p is sufficiently large, the D-parameters converge to values that are independent of p. In order to demonstrate this convergence, a more complex Mueller matrix was generated by multiplying an elliptical retarder of randomly generated input retardance vector by a partial depolarizer PD of the form
PD(a,b,c)=(10000a0000b0000c).
(32)
The randomly generated elliptical retarder with retardance vector {δH,δ45,δR} = {0.210,0.033,1.003} was multiplied by PD(0.1,0.2,0.3), generating the Mueller matrix
ED(0.210,0.033,1.003)PD(0.1,0.2,0.3)=(0.4940.0100.0030.4190.1040.0020.00.0290.0160.00.00.0050.4950.0100.0030.142).
(33)
Its Mueller matrix root decomposition parameters were then calculated for varying values of p. Figure 4 shows the convergence of the norm of the matrix roots retardance parameters (D4, D5, and D6). The norm of the matrix roots retardance parameters converges to a steady value of approximately 1.07 following the 104th root. The convergence of the norm of the matrix roots diagonal depolarization parameters (D13, D14, and D15) behaves in a strikingly similar manner, as shown in Fig. 3.

Fig. 3 The norm of the depolarizing retarder’s diagonal depolarization parameters (D13, D14, and D15) converges to a steady value of approximately 3.03 beyond the 104th root.
Fig. 4 The norm of the depolarizing retarder’s matrix roots retardance parameters (D4, D5, and D6) converges to a steady value of approximately 1.07 beyond the 104th root.

5.1.2. Numerical precision

Fig. 5 Relative error for Mathematica’s matrix root calculation.

Based on the behaviors documented in Figs. 2 through 5, a good choice for p is on the order of 105. This choice balances the relative error generated from the root calculation while achieving convergence of its matrix root polarization properties. Relative error of 5 · 10−12 is achieved with p = 105 for the randomly generated Mueller matrix MR.

6. Algorithm and flow chart

An algorithm to calculate the matrix roots decomposition incorporating the results from sections 4 and 5 follows.

Fig. 6 Matrix Roots algorithm flow chart.

7. Statistical algorithm implementation

The matrix roots decomposition algorithm from section 6 was implemented on 76,336 randomly generated, non-singular, physical Mueller matrices with no real negative eigenvalues. By definition, all of these Mueller matrices have principal roots. The statistical analysis resulting from this implementation yields information about the range in values of the depolarizing matrix roots parameters for a large quantity of random Mueller matrices.

The Mueller matrices are generated by a brute-force numerical method. A four-by-four matrix is generated by setting the m0,0 element to a value of one, and the other fifteen matrix elements are uniformly randomly distributed between negative one and one. If the matrix is nonphysical or has negative real eigenvalues, it is discarded. The Mueller matrix roots parameters were calculated for all of the remaining matrices. 76,336 physical, non-singular Mueller matrices with no real negative eigenvalues were found from the set of 109 randomly generated matrices.

Figures 7 through 9 show the distribution of the depolarizing matrix roots parameters. The distributions for the amplitude depolarization parameters (D7, D8, and D9) are entirely overlapping and largely lie within the range of −1 to 1, with a full-width half maximum (FWHM) of 0.6. The distributions for the phase depolarization parameters (D10, D11, and D12) also overlap entirely and range mostly between −2 and 2, with a FWHM of 0.8. The distributions of the diagonal depolarization parameters D13 and D14 overlap and have a very similar distribution to the phase depolarization parameters, with a range primarily between −2 and 2 and FWHM of 0.8. D15 has a distinct distribution. It has a hard limit at zero, as it cannot have a negative value, since f13, f14, and f15 are always positive. Its distribution cuts off near 4, with a FWHM of 1.1. Out of all the depolarizing matrix roots parameters, it has the only non-symmetric distribution.

Fig. 7 A histogram of the matrix root amplitude depolarization values for 76,336 randomly generated physical Mueller matrices.
Fig. 8 A histogram of the matrix root phase depolarization values for 76,336 randomly generated physical Mueller matrices.
Fig. 9 A histogram of the matrix root diagonal depolarization values for 76,336 randomly generated physical Mueller matrices.

8. Conclusion

Computational issues involved in applying the Mueller matrix roots decomposition have been addressed. The definition of the pth matrix root is reviewed, along with a brief discussion of the most common methods of calculating the pth principal matrix root. Our study indicates that the decomposition is optimized around p = 105 for a 53 digit binary double precision calculation, in consideration of numerical accuracy and noise as well as parameter convergence. Practical values for the roots of singular Mueller matrices can be obtained through perturbing them to nearby diattenuating matrices. Similarly, Mueller matrices with a half wave of retardance can be evaluated by perturbing their retardance from a half wave, without changing the retarder form. An algorithm is provided which incorporates the computational considerations involved in calculating the matrix roots decomposition. Finally, the algorithm is implemented to perform a statistical analysis on a large set of randomly generated Mueller matrices in order to yield insight on the typical ranges of the matrix roots parameters for physical Mueller matrices.

References and links

1.

S.-Y. Lu and R. A. Chipman, “Interpretation of mueller matrices based on polar decomposition,” J. Opt. Soc. Am. A 13, 1106–1113 (1996). [CrossRef]

2.

R. Ossikovski, “Analysis of depolarizing mueller matrices through a symmetric decomposition,” J. Opt. Soc. Am. A 26, 1109–1118 (2009). [CrossRef]

3.

S. R. Cloude, “Conditions for the physical realisability of matrix operators in polarimetry,” Proc. SPIE 1166, 177–185 (1989).

4.

R. A. Chipman, Handbook of Optics, Vol. 1 of Mueller Matrices (McGraw Hill, 2009), 3rd ed.

5.

R. C. Jones, “A new calculus for the treatment of optical systems. vii. properties of the n-matrices,” J. Opt. Soc. Am. 38, 671–683 (1948). [CrossRef]

6.

N. J. Higham, Functions of Matrices: Theory and Computation (Society for Industrial and Applied Mathematics, 2008). [CrossRef]

7.

D. A. Bini, N. J. Higham, and B. Meini, “Algorithms for the matrix pth root,” Num. Algor. 39, 349–378 (2005). [CrossRef]

8.

C.-H. Guo, “On newton’s method and halley’s method for the principal pth root of a matrix,” Linear Algebra Appl. 432, 1905–1922 (2010). [CrossRef]

9.

R. D. Skeel and J. B. Keiper, Elementary Numerical Computing with Mathematica, Chapter 2 (McGraw Hill, 1993).

10.

R. A. Chipman, “Depolarization index and the average degree of polarization,” Appl. Opt. 44, 2490–2495 (2005). [CrossRef] [PubMed]

11.

J. J. Gil and E. Bernabeu, “Depolarization and polarization indices of an optical system,” Opt. Acta 33, 185–189 (1986). [CrossRef]

12.

H. D. Noble, “Mueller matrix roots,” Doctoral dissertation (2011).

13.

N. J. Higham, “Stable iterations for the matrix square root,” Num. Algor. 15, 227–242 (1997). [CrossRef]

14.

M. H. Smith, “Optimization of a dual-rotating-retarder mueller matrix imaging polarimeter,” Appl. Opt. 41, 2488–2493 (2002). [CrossRef] [PubMed]

15.

S.-Y. Lu and R. A. Chipman, “Homogeneous and inhomogeneous jones matrices,” J. Opt. Soc. Am. A 11, 766–773 (1994). [CrossRef]

16.

D. Goldstein, Polarized Light, 2nd ed. (Marcel Dekker, 2003). [CrossRef]

17.

R. Barakat, “Conditions for the physical realizability of polarization matrices characterizing passive systems,” J. Mod. Opt. 34, 1535–1544 (1987). [CrossRef]

18.

K. M. Twietmeyer, R. A. Chipman, A. E. Elsner, Y. Zhao, and D. VanNasdale, “Mueller matrix retinal imager with optimized polarization conditions,” Opt. Express 16, 21339–21354 (2008). [CrossRef] [PubMed]

OCIS Codes
(120.5410) Instrumentation, measurement, and metrology : Polarimetry
(260.5430) Physical optics : Polarization
(290.5855) Scattering : Scattering, polarization
(240.2130) Optics at surfaces : Ellipsometry and polarimetry

ToC Category:
Instrumentation, Measurement, and Metrology

History
Original Manuscript: August 2, 2011
Revised Manuscript: October 13, 2011
Manuscript Accepted: October 21, 2011
Published: December 19, 2011

Citation
H. D. Noble and R. A. Chipman, "Mueller matrix roots algorithm and computational considerations," Opt. Express 20, 17-31 (2012)
http://www.opticsinfobase.org/oe/abstract.cfm?URI=oe-20-1-17


Sort:  Author  |  Year  |  Journal  |  Reset  

References

  1. S.-Y. Lu and R. A. Chipman, “Interpretation of mueller matrices based on polar decomposition,” J. Opt. Soc. Am. A13, 1106–1113 (1996). [CrossRef]
  2. R. Ossikovski, “Analysis of depolarizing mueller matrices through a symmetric decomposition,” J. Opt. Soc. Am. A26, 1109–1118 (2009). [CrossRef]
  3. S. R. Cloude, “Conditions for the physical realisability of matrix operators in polarimetry,” Proc. SPIE1166, 177–185 (1989).
  4. R. A. Chipman, Handbook of Optics, Vol. 1 of Mueller Matrices (McGraw Hill, 2009), 3rd ed.
  5. R. C. Jones, “A new calculus for the treatment of optical systems. vii. properties of the n-matrices,” J. Opt. Soc. Am.38, 671–683 (1948). [CrossRef]
  6. N. J. Higham, Functions of Matrices: Theory and Computation (Society for Industrial and Applied Mathematics, 2008). [CrossRef]
  7. D. A. Bini, N. J. Higham, and B. Meini, “Algorithms for the matrix pth root,” Num. Algor.39, 349–378 (2005). [CrossRef]
  8. C.-H. Guo, “On newton’s method and halley’s method for the principal pth root of a matrix,” Linear Algebra Appl.432, 1905–1922 (2010). [CrossRef]
  9. R. D. Skeel and J. B. Keiper, Elementary Numerical Computing with Mathematica, Chapter 2 (McGraw Hill, 1993).
  10. R. A. Chipman, “Depolarization index and the average degree of polarization,” Appl. Opt.44, 2490–2495 (2005). [CrossRef] [PubMed]
  11. J. J. Gil and E. Bernabeu, “Depolarization and polarization indices of an optical system,” Opt. Acta33, 185–189 (1986). [CrossRef]
  12. H. D. Noble, “Mueller matrix roots,” Doctoral dissertation (2011).
  13. N. J. Higham, “Stable iterations for the matrix square root,” Num. Algor.15, 227–242 (1997). [CrossRef]
  14. M. H. Smith, “Optimization of a dual-rotating-retarder mueller matrix imaging polarimeter,” Appl. Opt.41, 2488–2493 (2002). [CrossRef] [PubMed]
  15. S.-Y. Lu and R. A. Chipman, “Homogeneous and inhomogeneous jones matrices,” J. Opt. Soc. Am. A11, 766–773 (1994). [CrossRef]
  16. D. Goldstein, Polarized Light, 2nd ed. (Marcel Dekker, 2003). [CrossRef]
  17. R. Barakat, “Conditions for the physical realizability of polarization matrices characterizing passive systems,” J. Mod. Opt.34, 1535–1544 (1987). [CrossRef]
  18. K. M. Twietmeyer, R. A. Chipman, A. E. Elsner, Y. Zhao, and D. VanNasdale, “Mueller matrix retinal imager with optimized polarization conditions,” Opt. Express16, 21339–21354 (2008). [CrossRef] [PubMed]

Cited By

Alert me when this paper is cited

OSA is able to provide readers links to articles that cite this paper by participating in CrossRef's Cited-By Linking service. CrossRef includes content from more than 3000 publishers and societies. In addition to listing OSA journal articles that cite this paper, citing articles from other participating publishers will also be listed.


« Previous Article  |  Next Article »

OSA is a member of CrossRef.

CrossCheck Deposited