## 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 *p*th matrix root is optimized around *p* = 10^{5} 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]. 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]. The diattenuation, retardance, and depolarization parameters calculated by these methods depend on the arbitrary order of the decomposition. The Cloude additive decomposition [3] 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] 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] to Mueller matrices by dividing the Mueller matrix **M** into *p* infinitesimal slices:

**N**is the Mueller matrix for one small slice. When

*p*becomes very large, the polarization properties of $\sqrt[p]{\mathbf{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

*p*th root of a large class of Mueller matrices approaches the identity matrix

**I**, In this study, the class of Mueller matrices that obey Eq. (2) are called

*uniform*Mueller matrices. The Mueller matrix roots decomposition from [4] is applicable to this subset of Mueller matrices, with the exception of the non-uniform special cases highlighted in section 4.

The calculation of the principal *p*th root of square matrices is an extensively studied subject [6–8], and numerical accuracy and noise are well understood within the field of numerical computing [9]. 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] to be order-independent, the *p*th root must be optimized with all of these topics in mind. This paper reviews and updates the Mueller matrix roots decomposition from [4], 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 *p*th principal root of **M**, where *p* is some large integer (typically 10^{5}),

*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*_{0} through *d*_{15} are defined from the symmetric and antisymmetric parts of N:

*d*

_{13},

*d*

_{14}, and

*d*

_{15}are solved from the first-order generator products in terms of the parameters

*f*

_{13},

*f*

_{14}, and

*f*

_{15}: The infinitesimal polarization parameters

*d*

_{0}through

*d*

_{15}are rescaled by

*p*to produce the matrix roots parameters

*D*

_{0}through

*D*

_{15}:

*D*

_{0}through

*D*

_{15}parameterize the sixteen degrees of freedom of

**M**.

There are three matrix roots parameters for diattenuation, *D*_{1}, *D*_{2}, *D*_{3}, and three matrix roots parameters for retardance, *D*_{4}, *D*_{5}, *D*_{6}. 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*_{7} through *D*_{15}, describe depolarizing effects - one for each of the nine depolarizing degrees of freedom described in [4]. 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*_{7}, *D*_{8}, and *D*_{9}, 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]. They are named amplitude depolarization because they depolarize and affect the flux of an incident Stokes vector. The terms for phase depolarization, *D*_{10}, *D*_{11}, and *D*_{12}, 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*_{13}, *D*_{14}, and *D*_{15} are named diagonal depolarization because they lie on the matrix diagonal. *D*_{15} expresses the overall isotropic depolarizing power, since *D*_{15} reduces the degree of polarization of all incident Stokes vectors equally, independent of the Stokes vector’s location on the Poincare sphere. *D*_{13} expresses the relative strength between the diagonal depolarization on the two linear axes (horizontal/vertical and 45°/135°). *D*_{14} expresses the relative strength between linear and circular diagonal depolarization. The diagonal depolarization parameters have been modified from [4] so that only the diagonal depolarization parameter *D*_{15} changes the Gil-Bernabeu depolarization index [10,11], while *D*_{0} through *D*_{14} do not [12]. Table 1 lists the parameters *D*_{1} through *D*_{15} and categorizes them into their corresponding families and axes.

## 3. *p*th root

This section discusses the definition of the principal *p*th matrix root along with relevant examples and common methods of calculating the *p*th 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}*p*th roots [6]. So long as the matrix

**A**∈ ℂ

*has no negative, real eigenvalues, there is a unique*

^{n×n}*p*th root of

**A**whose eigenvalues’ arguments lie between −

*π*/

*p*and

*π*/

*p*, and that unique root is defined as the principal root of

**A**[6]. 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 *p*th root of the resulting upper triangular factor using various (stable) recursive formulae [13]. Newton’s method calculates the *p*th root of **A** using an iterative approach [6]. The Newton method is largely considered to have poor convergence and stability properties [14], and the Schur method from [14] is the numerically stable benchmark against which other methods are often compared [7]. The Schur-Newton algorithm applies iterative computations to the upper triangular matrices from the Schur decomposition [6]. Many of these algorithms are available in the MATLAB Matrix Computation Toolbox [7].

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,

*λ*are the eigenvalues of

_{i}**A**, the columns of

**Z**are its eigenvectors, and diag(

*λ*) is a diagonal matrix with its

_{i}*i*th diagonal element equal to

*λ*. Then the matrix root of

_{i}**A**is calculated as 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].

## 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] is always idempotent, the matrix roots decomposition will break down at the point where the *p*th root is calculated. For example, the formula for a homogeneous linear polarizer as a function of orientation *θ* is

A homogeneous polarizer has two orthogonal, physical eigenpolarizations with eigenvalues (*T _{max}*, 0), where

*T*is the maximum transmission, and the minimum transmission

_{max}*T*is 0. The polarizer can be replaced with a nearby uniform Mueller matrix by adjusting the maximum and minimum transmission by a small number

_{min}*ɛ*. 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*

_{45},

*d*} [1, 4], where

_{R}*η*and

*θ*are then plugged in to the general formula for an elliptical diattenuator [4, 16]. For an ideal polarizer,

*T*= 1 and

_{max}*T*= 0. The

_{min}*p*th 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,

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*p*th root of an ideal homogeneous linear retarder

**LR**(

*δ*,

*θ*) with retardance

*δ*at orientation

*θ*is 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°,

*λ*= {−1,−1,1,1}), so the half wave retarder has no principal

*p*th root. For any value of

*p*, two eigenvalues of

**HWR**

^{1/p}are (−1)

*with arguments of −*

^{p}*π*/

*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

*p*th 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},

*δ*} [1, 4] as follows,

_{R}*δ*, orientation

*θ*, and latitude

*η*can be written in terms of three linear retarders as

**LR**(

*δ*,

*θ*) with retardance

*δ*and an orientation of

*θ*is [4]

*θ*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:

It is also possible to calculate the root in Eq. (24) without perturbation of the retardance,

#### 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*

*λ*= {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] is first calculated as an intermediate step,

**M**

_{Δ}is a depolarizing Mueller matrix,

**M**

*is a pure retarder, and*

_{R}**M**

*is a pure diattenuator. The half-wave retarder found in the*

_{D}**M**

*term is perturbed to the nearby retarder*

_{R}**M**′

*with (*

_{R}*π*–

*ɛ*) retardance as shown in equation 24. Then the terms are recombined to form

**M**′, where The eigenvalues of

**M**′ are no longer negative, since

**M**

_{Δ}and

**M**

*always have positive eigenvalues [1], and the eigenvalues of*

_{D}**M**′

*are positive. Therefore the uniform*

_{R}*p*th 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] is first calculated as an intermediate step,

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**′

*is substituted into Eq. (29) in place of*

_{D}**M**

*, and the three matrices are recombined to form a modified, non-singular matrix*

_{D}**M**′: Then the

*p*th 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 *p*th root of a retarder Mueller matrix with retardance *δ* results in a *p*th principal matrix root with retardance *δ*/*p*. Therefore a Mueller matrix of an ideal elliptical retarder with retarder vector {*δ _{H}*,

*δ*

_{45},

*δ*} should have matrix roots parameters of

_{R}*D*

_{4}=

*δ*,

_{H}*D*

_{5}=

*δ*

_{45}, and

*D*

_{6}=

*δ*. In order to evaluate the correct choice of

_{R}*p*in consideration of this criteria, the matrix roots retardance parameters (

*D*

_{4}through

*D*

_{6}) were calculated (using Mathematica’s default double-precision machine arithmetic) for a pure elliptical retarder with retardance parameters {

*δ*,

_{H}*δ*

_{45},

*δ*} = {0.294,0.302,0.997} for different values of root order

_{R}*p*. The relative error Δ

*was then calculated according to the following expression,*

_{x}Figure 2 shows the relative error between the input retardance vector {*δ _{H}*,

*δ*

_{45},

*δ*} and the corresponding matrix roots retardance vector

_{R}*D*

_{4}through

*D*

_{6}for different choices of

*p*for the

*p*th root. In Fig. 2, the relative retardance error converges to a minimum value just beyond the 10

^{5}th root. For large

*p*, the relative error increases, due to numerical rounding and the loss of precision associated with machine arithmetic.

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

*δ*,

_{H}*δ*

_{45},

*δ*} = {0.210,0.033,1.003} was multiplied by

_{R}**PD**(0.1,0.2,0.3), generating the Mueller matrix

*p*. Figure 4 shows the convergence of the norm of the matrix roots retardance parameters (

*D*

_{4},

*D*

_{5}, and

*D*

_{6}). The norm of the matrix roots retardance parameters converges to a steady value of approximately 1.07 following the 10

^{4}th root. The convergence of the norm of the matrix roots diagonal depolarization parameters (

*D*

_{13},

*D*

_{14}, and

*D*

_{15}) behaves in a strikingly similar manner, as shown in Fig. 3.

### 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]. 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].

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

*, for*

^{k}*k*between 1 and 14. This relative error was calculated as follows,

Based on the behaviors documented in Figs. 2 through 5, a good choice for *p* is on the order of 10^{5}. 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* = 10^{5} 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], if the matrix is nonphysical, it is recommended to use a force-physical routine [17,18] 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^{5}th matrix root of the Mueller matrix. (This choice of *p* = 10^{5} 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*_{0} through *d*_{15} are determined from the principal root according to Eq. (4). Step 8 rescales the infinitesimal roots parameters by *p* = 10^{5}, resulting in the desired matrix roots parameters *D*_{0} through *D*_{15}.

## 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^{9} randomly generated matrices.

Figures 7 through 9 show the distribution of the depolarizing matrix roots parameters. The distributions for the amplitude depolarization parameters (*D*_{7}, *D*_{8}, and *D*_{9}) 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*_{10}, *D*_{11}, and *D*_{12}) also overlap entirely and range mostly between −2 and 2, with a FWHM of 0.8. The distributions of the diagonal depolarization parameters *D*_{13} and *D*_{14} 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*_{15} has a distinct distribution. It has a hard limit at zero, as it cannot have a negative value, since *f*_{13}, *f*_{14}, and *f*_{15} 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.

## 8. Conclusion

Computational issues involved in applying the Mueller matrix roots decomposition have been addressed. The definition of the *p*th matrix root is reviewed, along with a brief discussion of the most common methods of calculating the *p*th principal matrix root. Our study indicates that the decomposition is optimized around *p* = 10^{5} 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]