Mueller matrix roots algorithm and computational considerations

doi:10.1364/OE.20.000017

« 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)

Journal Search
Article Lookup

Article Tools

Citations

Alert me when this article is cited
Export Citation/Save

Full Text
Article Info
References (18)
Cited By
Figures (9)
Metrics

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 = 10⁵ 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.

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

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

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

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

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

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 = \sqrt[p]{M},

(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

\sqrt[p]{M}

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,

lim_{p \to \infty} \sqrt[p]{M} = 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

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

N. J. Higham, Functions of Matrices: Theory and Computation (Society for Industrial and Applied Mathematics, 2008). [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]

], and numerical accuracy and noise are well understood within the field of numerical computing [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

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

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 10⁵),

N = \sqrt[p]{M} .

(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 d₀ through d₁₅ are defined from the symmetric and antisymmetric parts of N:

N = d_{0} (\begin{matrix} 1 & d_{1} + d_{7} & d_{2} + d_{8} & d_{3} + d_{9} \\ d_{1} - d_{7} & 1 - f_{13} & d_{6} + d_{12} & - d_{5} + d_{11} \\ d_{2} - d_{8} & - d_{6} + d_{12} & 1 - f_{14} & d_{4} + d_{10} \\ d_{3} - d_{9} & d_{5} + d_{11} & - d_{4} + d_{10} & 1 - f_{15} \end{matrix}),

(4)

where d₁₃, d₁₄, and d₁₅ are solved from the first-order generator products in terms of the parameters f₁₃, f₁₄, and f₁₅:

d_{13} = \frac{f_{14} - f_{13}}{2}

(5)

d_{14} = \frac{2}{\sqrt{3}} (- 2 f_{15} + f_{14} + f_{13})

(6)

d_{15} = \frac{1}{\sqrt{6}} (f_{15} + f_{14} + f_{13}) .

(7)

The infinitesimal polarization parameters d₀ through d₁₅ are rescaled by p to produce the matrix roots parameters D₀ through D₁₅:

D_{n} = p \cdot d_{n}, n = 0, 1, 2, \dots, 15 .

(8)

D₀ through D₁₅ parameterize the sixteen degrees of freedom of M.

There are three matrix roots parameters for diattenuation, D₁, D₂, D₃, and three matrix roots parameters for retardance, D₄, D₅, D₆. 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, D₇ through D₁₅, describe depolarizing effects - one for each of the nine depolarizing degrees of freedom described in [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, D₇, D₈, and D₉, 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

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, D₁₀, D₁₁, and D₁₂, 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. D₁₃, D₁₄, and D₁₅ are named diagonal depolarization because they lie on the matrix diagonal. D₁₅ expresses the overall isotropic depolarizing power, since D₁₅ reduces the degree of polarization of all incident Stokes vectors equally, independent of the Stokes vector’s location on the Poincare sphere. D₁₃ expresses the relative strength between the diagonal depolarization on the two linear axes (horizontal/vertical and 45°/135°). D₁₄ expresses the relative strength between linear and circular diagonal depolarization. The diagonal depolarization parameters have been modified from [4

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

] so that only the diagonal depolarization parameter D₁₅ changes the Gil-Bernabeu depolarization index [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]

], while D₀ through D₁₄ do not [12

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

]. Table 1 lists the parameters D₁ through D₁₅ 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.

Property	Horizontal/Vertical	45°/135°	Right/Left Circular
Diattenuation	D₁	D₂	D₃
Retardance	D₄	D₅	D₆
Amplitude Depolarization	D₇	D₈	D₉
Phase Depolarization	D₁₀	D₁₁	D₁₂
	Relative Linear	Relative Linear or Circular	Isotropic
Diagonal Depolarization	D₁₃	D₁₄	D₁₅

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 p^s pth roots [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

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 A^1/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

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

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

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

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

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

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

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 = Z diag (λ_{i}) Z^{- 1},

(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

A^{1 / p} = Z diag {(λ_{i})}^{1 / p} Z^{- 1} .

(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 A^1/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

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

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 (θ) = \frac{1}{2} (\begin{matrix} 1 & cos 2 θ & sin 2 θ & 0 \\ cos 2 θ & {cos}^{2} 2 θ & \frac{1}{2} sin 4 θ & 0 \\ sin 2 θ & \frac{1}{2} sin 4 θ & {sin}^{2} 2 θ & 0 \\ 0 & 0 & 0 & 0 \end{matrix}) .

(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 (T_max, 0), where T_max is the maximum transmission, and the minimum transmission T_min 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 (T_max – ɛ,ɛ).

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 {d_H, d₄₅, d_R} [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

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

], where

θ = \frac{1}{2} arctan \frac{d_{H}}{d_{45}}, and

(14)

η = arcsin \frac{d_{R}}{\sqrt{d_{H}^{2} + d_{45}^{2} + d_{R}^{2}}} .

(15)

η and θ are then plugged in to the general formula for an elliptical diattenuator [4

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

, 16

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

]. For an ideal polarizer, T_max = 1 and T_min = 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 °) = (\begin{array}{l} 1 & 0 & 0 & 0 \\ 0 & 1 & 0 & 0 \\ 0 & 0 & - 1 & 0 \\ 0 & 0 & 0 & - 1 \end{array}),

(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 HWR^1/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, δ₄₅, δ_R} [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

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

] as follows,

θ = \frac{1}{2} arctan \frac{δ_{H}}{δ_{45}}, and

(19)

η = arcsin \frac{δ_{R}}{\sqrt{δ_{H}^{2} + δ_{45}^{2} + δ_{R}^{2}}} .

(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

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

]

LR (δ, θ) = (\begin{matrix} 1 & 0 & 0 & 0 \\ 0 & {cos}^{2} 2 θ + cos δ {sin}^{2} 2 θ & {sin}^{2} \frac{δ}{2} sin 4 θ & - sin δ sin 2 θ \\ 0 & {sin}^{2} \frac{δ}{2} sin 4 θ & cos δ {cos}^{2} 2 θ + {sin}^{2} 2 θ & cos 2 θ sin δ \\ 0 & sin δ sin 2 θ & - cos 2 θ sin δ & cos δ \end{matrix}) .

(22)

After θ and η are calculated from the retardance vector, the half-wave retarder is perturbed by some small number (ɛ ≈ 10⁻⁷) 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 {(η, θ, δ = π)}^{\frac{1}{p}} = LR (η, θ + π / 4) LR {((π - ɛ), θ)}^{\frac{1}{p}} LR (η, θ - π / 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 {(η, θ)}^{\frac{1}{p}} = LR (ɛ, θ + π / 4) LR (\frac{π}{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

(\begin{matrix} 1 & 0 & 0 & 0 \\ 0 & 1 & 0 & 0 \\ 0 & 0 & - 1 & 0 \\ 0 & 0 & 0 & - 1 \end{matrix}) (\begin{matrix} 1 & 0 & 0 & 0 \\ 0 & a & 0 & 0 \\ 0 & 0 & b & 0 \\ 0 & 0 & 0 & c \end{matrix}) = (\begin{matrix} 1 & 0 & 0 & 0 \\ 0 & a & 0 & 0 \\ 0 & 0 & - b & 0 \\ 0 & 0 & 0 & - c \end{matrix}),

(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 b ≠ c, 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

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_{Δ} M_{R} M_{D} .

(27)

M_Δ is a depolarizing Mueller matrix, M_R is a pure retarder, and M_D is a pure diattenuator. The half-wave retarder found in the M_R term is perturbed to the nearby retarder M′_R with (π–ɛ) retardance as shown in equation 24. Then the terms are recombined to form M′, where

M^{'} = M_{Δ} M_{R}^{'} M_{D} .

(28)

The eigenvalues of M′ are no longer negative, since M_Δ and M_D always have positive eigenvalues [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 M′_R 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

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_{Δ} M_{R} M_{D} .

(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 M′_D is substituted into Eq. (29) in place of M_D, and the three matrices are recombined to form a modified, non-singular matrix M′:

M^{'} = M_{Δ} M_{R} M_{D}^{'} .

(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,δ₄₅,δ_R} should have matrix roots parameters of D₄ = δ_H, D₅ = δ₄₅, and D₆ = δ_R. In order to evaluate the correct choice of p in consideration of this criteria, the matrix roots retardance parameters (D₄ through D₆) were calculated (using Mathematica’s default double-precision machine arithmetic) for a pure elliptical retarder with retardance parameters {δ_H,δ₄₅,δ_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} = \frac{\sqrt{{(δ_{H} - D_{4})}^{2} + {(δ_{45} - D_{5})}^{2} + {(δ_{R} - D_{6})}^{2}}}{δ} .

(31)

Figure 2 shows the relative error between the input retardance vector {δ_H,δ₄₅,δ_R} and the corresponding matrix roots retardance vector D₄ through D₆ for different choices of p for the pth root. In Fig. 2, the relative retardance error converges to a minimum value just beyond the 10⁵th 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 10⁵th 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) = (\begin{matrix} 1 & 0 & 0 & 0 \\ 0 & a & 0 & 0 \\ 0 & 0 & b & 0 \\ 0 & 0 & 0 & c \end{matrix}) .

(32)

The randomly generated elliptical retarder with retardance vector {δ_H,δ₄₅,δ_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) \cdot PD (0.1, 0.2, 0.3) = (\begin{matrix} 0.494 & 0.010 & 0.003 & 0.419 \\ 0.104 & 0.002 & 0.0 & 0.029 \\ 0.016 & 0.0 & 0.0 & 0.005 \\ 0.495 & 0.010 & 0.003 & 0.142 \end{matrix}) .

(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 (D₄, D₅, and D₆). The norm of the matrix roots retardance parameters converges to a steady value of approximately 1.07 following the 10⁴th root. The convergence of the norm of the matrix roots diagonal depolarization parameters (D₁₃, D₁₄, and D₁₅) behaves in a strikingly similar manner, as shown in Fig. 3.

Fig. 3 The norm of the depolarizing retarder’s diagonal depolarization parameters (D₁₃, D₁₄, and D₁₅) converges to a steady value of approximately 3.03 beyond the 10⁴th root.

Fig. 4 The norm of the depolarizing retarder’s matrix roots retardance parameters (D₄, D₅, and D₆) converges to a steady value of approximately 1.07 beyond the 10⁴th root.

5.1.2. Numerical precision

This section addresses issues related to the numerical precision of the matrix root calculation, using Mathematica’s built-in commands. Computational software programs such as Mathematica and Matlab use floating-point numbers with machine precision by default [9

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

]. The value of machine precision that produced the results included here is 15.96 digits, which corresponds to a 53 digit binary double precision number with a mantissa [9

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

Figure 5 shows the relative error (on a log scale) of the matrix root calculation on a randomly generated Mueller matrix M_R as a function of root order p = 10^k, for k between 1 and 14. This relative error was calculated as follows,

Δ_{root} = \frac{‖ M_{R} - {(\sqrt[p]{M_{R}})}^{p} ‖}{2 ‖ M_{R} ‖} .

(34)

The relative error increases linearly with each operation, so it is recommended to balance this with the optimization of the convergence properties discussed in section 5.1.1. The linear increase in error with each numerical operation is independent of the choice of Mueller matrix.

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 10⁵. 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⁻¹² is achieved with p = 10⁵ for the randomly generated Mueller matrix M_R.

6. Algorithm and flow chart

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

Figure 6 shows a flow chart for the algorithm, beginning with an input Mueller matrix. If the Mueller matrix is not physical, the decomposition parameters will not be meaningful. Since many measured Mueller matrices are slightly nonphysical [16

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

], if the matrix is nonphysical, it is recommended to use a force-physical routine [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]

] in Step 2. Step 3 tests for the special cases discussed in section 4 - Mueller matrices with a half wave of retardance and singular Mueller matrices. If the Mueller matrix tests positive it is perturbed according to the procedure in section 4, and the algorithm resumes. Step 4 tests for negative real eigenvalues, and if present, the decomposition is aborted. Step 5 calculates the 10⁵th matrix root of the Mueller matrix. (This choice of p = 10⁵ was discussed in section 5.) The matrix root can be calculated using a builtin matrix root algorithm (such as Mathematica’s matrix power function), or with any principal matrix root algorithm such as those discussed in section 3.1. The principal root algorithms discussed in section 3.1 can fail to converge to a solution, particularly for such a high root order. As discussed in section 3.1, for any matrix with imaginary eigenvalues, computational rounding errors can lead to spurious imaginary parts. These may be discarded so long as they are small enough not to affect the accuracy of the calculation. Step 6 discards these spurious imaginary parts. In step 7, the parameters d₀ through d₁₅ are determined from the principal root according to Eq. (4). Step 8 rescales the infinitesimal roots parameters by p = 10⁵, resulting in the desired matrix roots parameters D₀ through D₁₅.

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 m_0,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 10⁹ randomly generated matrices.

Figures 7 through 9 show the distribution of the depolarizing matrix roots parameters. The distributions for the amplitude depolarization parameters (D₇, D₈, and D₉) 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 (D₁₀, D₁₁, and D₁₂) also overlap entirely and range mostly between −2 and 2, with a FWHM of 0.8. The distributions of the diagonal depolarization parameters D₁₃ and D₁₄ 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. D₁₅ has a distinct distribution. It has a hard limit at zero, as it cannot have a negative value, since f₁₃, f₁₄, and f₁₅ 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 = 10⁵ 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

S.-Y. Lu and R. A. Chipman, “Interpretation of mueller matrices based on polar decomposition,” J. Opt. Soc. Am. A13, 1106–1113 (1996). [CrossRef]
R. Ossikovski, “Analysis of depolarizing mueller matrices through a symmetric decomposition,” J. Opt. Soc. Am. A26, 1109–1118 (2009). [CrossRef]
S. R. Cloude, “Conditions for the physical realisability of matrix operators in polarimetry,” Proc. SPIE1166, 177–185 (1989).
R. A. Chipman, Handbook of Optics, Vol. 1 of Mueller Matrices (McGraw Hill, 2009), 3rd ed.
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]
N. J. Higham, Functions of Matrices: Theory and Computation (Society for Industrial and Applied Mathematics, 2008). [CrossRef]
D. A. Bini, N. J. Higham, and B. Meini, “Algorithms for the matrix pth root,” Num. Algor.39, 349–378 (2005). [CrossRef]
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]
R. D. Skeel and J. B. Keiper, Elementary Numerical Computing with Mathematica, Chapter 2 (McGraw Hill, 1993).
R. A. Chipman, “Depolarization index and the average degree of polarization,” Appl. Opt.44, 2490–2495 (2005). [CrossRef] [PubMed]
J. J. Gil and E. Bernabeu, “Depolarization and polarization indices of an optical system,” Opt. Acta33, 185–189 (1986). [CrossRef]
H. D. Noble, “Mueller matrix roots,” Doctoral dissertation (2011).
N. J. Higham, “Stable iterations for the matrix square root,” Num. Algor.15, 227–242 (1997). [CrossRef]
M. H. Smith, “Optimization of a dual-rotating-retarder mueller matrix imaging polarimeter,” Appl. Opt.41, 2488–2493 (2002). [CrossRef] [PubMed]
S.-Y. Lu and R. A. Chipman, “Homogeneous and inhomogeneous jones matrices,” J. Opt. Soc. Am. A11, 766–773 (1994). [CrossRef]
D. Goldstein, Polarized Light, 2nd ed. (Marcel Dekker, 2003). [CrossRef]
R. Barakat, “Conditions for the physical realizability of polarization matrices characterizing passive systems,” J. Mod. Opt.34, 1535–1544 (1987). [CrossRef]
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.

Click here to see a list of articles that cite this paper