Fact-checked by Grok 2 weeks ago

Polar decomposition

In linear algebra, the polar decomposition is a fundamental theorem that factorizes any complex matrix A \in \mathbb{C}^{m \times n} as A = UH, where U is a partial isometry (unitary in the case of square invertible matrices) and H is a positive semi-definite Hermitian matrix uniquely determined as the square root of A^H A.^[1] This decomposition separates the "rotational" or isometric component U from the "stretching" or positive scaling component H, providing insight into the geometric action of the linear transformation represented by A.^[2] For real matrices, the decomposition takes the form A = QP, with Q orthogonal and P positive semi-definite symmetric.^[3] The theorem guarantees the existence and uniqueness of this factorization under appropriate conditions, such as when the rank of A equals the dimension of the domain for uniqueness of U.^[4] It is closely related to the singular value decomposition (SVD), from which the polar decomposition can be derived by setting U = P Q^H and H = Q \Sigma Q^H for A = P \Sigma Q^H, where \Sigma is diagonal with non-negative entries.^[4] Originally introduced by Léon Autonne in 1902 for normal matrices and later generalized, the polar decomposition has roots in early 20th-century work on matrix canonical forms and has been extensively developed in operator theory for Hilbert spaces.^[4] Beyond pure mathematics, the polar decomposition finds applications in numerical linear algebra for orthogonalization and matrix sign functions,^[4] in quantum mechanics for quantum algorithms and measurement theory,^[5] and in continuum mechanics for decomposing deformation tensors into stretch and rotation parts.^[6] Computationally, it can be obtained via iterative methods like the Newton-Schulz iteration or directly from the SVD, with numerical stability ensured in finite-precision arithmetic.^[4]

Core Concepts

Definition

In linear algebra, for square matrices, the polar decomposition factorizes a matrix into a product involving a unitary matrix and a positive semi-definite matrix, analogous to separating rotational and scaling components of a linear transformation.^[7] For any square matrix A \in \mathbb{C}^{n \times n}, the right polar decomposition expresses A as A = U P, where U is a unitary matrix satisfying U^* U = I_n (with U^* denoting the conjugate transpose and I_n the identity matrix), and P is a positive semi-definite Hermitian matrix.^[8] Here, positive semi-definite means P = P^* (Hermitian) with all eigenvalues non-negative, and P is explicitly the unique positive semi-definite square root of A^* A, denoted P = (A^* A)^{1/2}, satisfying P^2 = A^* A.^[8] The left polar decomposition takes the form A = P' U', where P' is the unique positive semi-definite Hermitian square root of A A^*, so P' = (A A^*)^{1/2} with (P')^2 = A A^*, and U' is unitary.^[9] If A is invertible, both the right and left decompositions are unique.^[8]

Geometric Interpretation

The polar decomposition of a linear transformation represented by a matrix A \in \mathbb{R}^{n \times n} separates it into a product A = UP, where U is an orthogonal matrix encoding a rotation or reflection, and P is a positive semi-definite symmetric matrix representing a stretching or compression along principal axes.^[3] In the two-dimensional case, for a $2 \times 2 matrix, this decomposition intuitively captures how A distorts shapes by first applying an elliptical stretch via P and then a rigid rotation or reflection via U.^[10] A key visualization arises when considering the action on the unit circle in \mathbb{R}^2. The positive part P transforms the unit circle into an ellipse, with the semi-major and semi-minor axes aligned to the eigenvectors of P and lengths given by its eigenvalues (the principal stretches). The unitary part U then rotates or reflects this ellipse to its final position, matching the image of the unit circle under A. This process highlights the decomposition's role in isolating distortive and rigid components of the transformation.^[10] For a concrete example, consider the matrix A = \begin{pmatrix} 2 & 0 \\ 0 & 1 \end{pmatrix}, which stretches vectors along the x-axis by a factor of 2 while leaving the y-axis unchanged. Its polar decomposition is A = I \cdot P, where U = I (the identity matrix, implying no rotation) and P = A itself, as P is already positive semi-definite. Geometrically, this simply elongates the unit circle into an ellipse with axes of lengths 2 and 1, without further reorientation.^[2] In three dimensions, the interpretation extends analogously: U performs an orthogonal transformation, such as a rotation around an axis or a reflection, preserving volumes and angles, while P applies a positive semi-definite scaling that deforms a unit sphere into an ellipsoid (degenerate if singular) along its principal axes.^[6] The unitary component U acts as an isometry, maintaining lengths of vectors, whereas P encodes the singular values of A as the stretch factors along orthogonal directions, quantifying the transformation's distortive effect independent of orientation.^[3]

Algebraic Properties

General Properties

In the polar decomposition A = U P of a square complex matrix A, where U is unitary and P is positive semi-definite Hermitian, the factor U is isometric and thus has operator norm \|U\| = 1. Since \|A x\| = \|U (P x)\| = \|P x\| for all x, it follows that the operator norm satisfies \|A\| = \|P\|.^[1] The eigenvalues of P coincide with the singular values of A, as P = \sqrt{A^* A} and the singular values are the square roots of the eigenvalues of A^* A. The eigenvalues of the unitary factor U lie on the unit circle in the complex plane. The Frobenius norm of A relates directly to P via \|A\|_F^2 = \operatorname{trace}(A^* A) = \operatorname{trace}(P^2), since A^* A = P^2. A matrix A is invertible if and only if its positive factor P is invertible, as U is always invertible; in this case, the inverse is given by A^{-1} = P^{-1} U^*. If A = U P and B = V Q are polar decompositions, the polar decomposition of the product AB does not take the simple form U V \cdot P Q, but instead has positive factor \sqrt{B^* P^2 B}.

Relation to SVD

The singular value decomposition (SVD) provides an explicit means to construct the polar decomposition of a matrix A \in \mathbb{C}^{m \times n}, expressed as A = W \Sigma V^*, where W \in \mathbb{C}^{m \times m} and V \in \mathbb{C}^{n \times n} are unitary matrices, and \Sigma \in \mathbb{R}^{m \times n} is a rectangular diagonal matrix with non-negative real singular values \sigma_1 \geq \sigma_2 \geq \cdots \geq \sigma_r > 0 on the diagonal (with r = \operatorname{rank}(A)) and zeros elsewhere.^[11] From this SVD, the right polar decomposition A = U P follows by setting U = W V^* (a partial isometry with U^* U being the orthogonal projection onto the range of V) and P = V \Sigma V^* (a positive semi-definite Hermitian matrix), yielding U P = W V^* V \Sigma V^* = W \Sigma V^* = A. The singular values of A are precisely the eigenvalues of P. Similarly, the left polar decomposition A = P' U' is obtained with P' = W \Sigma W^* (positive semi-definite Hermitian) and U' = W V^* (partial isometry), satisfying P' U' = W \Sigma W^* W V^* = W \Sigma V^* = A.^[11]^[12] For real matrices A \in \mathbb{R}^{m \times n}, the SVD takes the form A = W \Sigma V^T with orthogonal W and V, leading to analogous real polar forms where U and U' are orthogonal (or partial isometries) and P, P' are symmetric positive semi-definite; the inherent sign ambiguities in the columns of W and V allow choices that align phases to ensure properties like positive determinant in the orthogonal factors when desired. Computationally, these relations imply that the polar decomposition arises as a direct byproduct of the SVD via simple matrix multiplications, whereas extracting the SVD from a known polar form necessitates the eigendecomposition of the positive semi-definite factor, which is comparably intensive to SVD computation itself.

Relation to Normal Matrices

A normal matrix A satisfies A A^* = A^* A, where A^* denotes the adjoint (conjugate transpose) of A. In the polar decomposition A = U P, with U unitary and P positive semidefinite, A is normal if and only if U and P commute, i.e., U P = P U.^[13] This commutativity implies that there exists a unitary matrix Q that simultaneously diagonalizes both U and P, so A is unitarily similar to a diagonal matrix via the same basis.^[13] Since P^2 = A^* A is positive semidefinite (hence normal) and A is normal, the polar factors share a common eigenbasis. By the spectral theorem for normal matrices, A = Q \Lambda Q^*, where Q is unitary and \Lambda = \operatorname{diag}(\lambda_1, \dots, \lambda_n) contains the eigenvalues of A. The positive part is then P = Q |\Lambda| Q^*, with |\Lambda| = \operatorname{diag}(|\lambda_1|, \dots, |\lambda_n|), and the unitary part is U = Q (\Lambda / |\Lambda|) Q^*, where \Lambda / |\Lambda| = \operatorname{diag}(e^{i \arg(\lambda_1)}, \dots, e^{i \arg(\lambda_n)}). For example, consider the normal matrix with eigenvalues $2e^{i\pi/4} and $3e^{i\pi/3}; then P has diagonal entries 2 and 3, while U has e^{i\pi/4} and e^{i\pi/3} in the shared eigenbasis. This structure aligns the polar decomposition directly with the spectral theorem, as the eigenvalues of A determine both the singular values (of P) and the phases (of U).^[13]

Existence and Construction

For Normal Operators

For normal operators on a Hilbert space, the spectral theorem facilitates both the existence and explicit construction of the polar decomposition. Specifically, a bounded normal operator A on a Hilbert space \mathcal{H} admits a spectral decomposition given by

A = \int_{\sigma(A)} \lambda \, dE(\lambda),

where \sigma(A) denotes the spectrum of A and E is the unique spectral measure (resolution of the identity) such that E(\cdot) is a projection-valued measure on the Borel subsets of \mathbb{C} with support in \sigma(A).^[14] The positive part P = |A| of the decomposition is constructed via functional calculus as

P = \int_{\sigma(A)} |\lambda| \, dE(\lambda).

To see that P serves as the modulus, note that the self-adjoint operator A^* A has spectral decomposition

A^* A = \int_{\sigma(A)} |\lambda|^2 \, dE(\lambda),

and since the function f(t) = \sqrt{t} (for t \geq 0) is continuous and positive, the functional calculus yields P as the unique positive square root of A^* A, ensuring P is positive and self-adjoint.^[14] The partial isometry U is then defined by

U = \int_{\sigma(A) \setminus \{0\}} \frac{\lambda}{|\lambda|} \, dE(\lambda)

on the closure of the range of P, with U extended by zero on the kernel of P (which corresponds to the support of E(\{0\})). This ensures U^* U = I on the range of P and U U^* = I on the range of A, verifying that U is a partial isometry with initial space \overline{\operatorname{ran} P} and final space \overline{\operatorname{ran} A}. The product A = U P follows directly from the multiplicativity of the functional calculus, as

U P = \left( \int_{\sigma(A) \setminus \{0\}} \frac{\lambda}{|\lambda|} \, dE(\lambda) \right) \left( \int_{\sigma(A)} |\lambda| \, dE(\lambda) \right) = \int_{\sigma(A)} \lambda \, dE(\lambda) = A

on the range, with both sides vanishing on the kernel.^[14] The uniqueness of this polar decomposition for normal operators stems from the uniqueness of the spectral measure E in the spectral theorem, which determines both P and U uniquely via the above integrals. Any other decomposition A = V Q with Q positive and V a partial isometry satisfying the range conditions would yield the same Q = |A| (as the unique positive square root) and thus the same V = A Q^{-1} (where defined), aligning with the functional calculus construction.^[14] In the finite-dimensional case, where A is a normal matrix on \mathbb{C}^n, the spectral theorem reduces to unitary diagonalization A = W D W^* with D diagonal containing the eigenvalues \lambda_i. The polar decomposition then aligns with P = W \operatorname{diag}(|\lambda_1|, \dots, |\lambda_n|) W^* and U = W \operatorname{diag}(e^{i \arg \lambda_1}, \dots, e^{i \arg \lambda_k}, 0, \dots, 0) W^*, where the first k entries correspond to the non-zero eigenvalues.

For Invertible Matrices

For an invertible matrix A \in \mathbb{C}^{n \times n}, the adjoint A^* A is positive definite Hermitian, ensuring that the principal matrix logarithm \log(A^* A) exists and is Hermitian, as the eigenvalues of A^* A are positive real numbers. This setup allows a constructive approach to the polar decomposition via the matrix exponential, which maps Hermitian matrices to positive definite ones. The positive definite factor P is constructed as P = \exp\left(\frac{1}{2} \log(A^* A)\right), which coincides with the unique positive definite square root (A^* A)^{1/2}. The unitary factor is then U = A P^{-1}. This yields the right polar decomposition A = U P, where P is positive definite and U is unitary. To verify unitarity, note that U^* = P^{-1} A^*, so U^* U = P^{-1} A^* A P^{-1}. Substituting A^* A = P^2 gives U^* U = P^{-1} P^2 P^{-1} = I. Similarly, U U^* = I follows by symmetry or direct computation. In the real case, for invertible A \in \mathbb{R}^{n \times n}, A^T A is positive definite symmetric, so \log(A^T A) is real symmetric and P = \exp\left(\frac{1}{2} \log(A^T A)\right) is positive definite symmetric, yielding an orthogonal U = A P^{-1} with U^T U = I. The construction ensures orthogonality. The sign of \det U is determined by \det A, as \det P > 0. An alternative construction in the real case parametrizes the orthogonal factor U via the Cayley transform, where skew-symmetric matrices are mapped to orthogonal ones without eigenvalue -1, providing a rational approximation useful for avoiding logarithmic computations in certain numerical contexts.

General Proof

One approach to proving the existence of the polar decomposition for an arbitrary square complex matrix A \in \mathbb{C}^{n \times n}, possibly singular, relies on continuity arguments from the invertible case. Approximate A by the invertible matrix A_\epsilon = A + \epsilon I for \epsilon > 0, which admits a polar decomposition A_\epsilon = U_\epsilon P_\epsilon where U_\epsilon is unitary and P_\epsilon is positive definite. As \epsilon \to 0^+, the sequence \{P_\epsilon\} converges in the operator norm to a positive semi-definite matrix P = \lim_{\epsilon \to 0^+} P_\epsilon = (A^* A)^{1/2}, the unique positive semi-definite square root of A^* A. Similarly, \{U_\epsilon\} converges to a partial isometry U such that U is unitary on the range of A (equivalently, the range of P) and U P = A, with U acting as zero on the kernel of A. This limit yields the polar decomposition A = U P, as the approximation ensures the factors remain bounded and the decomposition is continuous in the space of matrices. An alternative proof constructs the decomposition directly using the singular value decomposition (SVD) of A, which exists for any square matrix. The SVD states that A = V \Sigma W^* where V and W are unitary matrices and \Sigma is diagonal with non-negative real entries (the singular values of A). Define P = W \Sigma W^*, which is positive semi-definite, and U = V W^*, which is unitary. Then A = U P, establishing the polar decomposition without relying on limits. This construction highlights that the proof is non-constructive in the sense that it invokes the existence of the SVD, whose proof itself requires spectral theory for Hermitian matrices.^[12] The polar decomposition A = U P is unique in the following sense: the positive semi-definite factor P is uniquely determined as the positive square root of A^* A. The unitary factor U is unique on the range of P (or equivalently, the range of A), where it acts as an isometry mapping the range of P onto the range of A; on the orthogonal complement of the range of P, U may be chosen arbitrarily as long as it remains unitary overall.^[1] The construction preserves rank, with \operatorname{rank}(A) = \operatorname{rank}(P), since the non-zero singular values of A correspond exactly to the positive eigenvalues of P, and the kernel dimensions match via the partial isometry property of U.^[12]

Generalizations

Bounded Operators on Hilbert Spaces

In the context of bounded linear operators on a Hilbert space H, the polar decomposition generalizes the finite-dimensional matrix case by factoring a bounded operator A \in B(H) as A = U |A|, where |A| = \sqrt{A^* A} is the unique positive square root of the positive self-adjoint operator A^* A, and U is a partial isometry with initial projection onto the closure of the range of |A|.^[15]^[16] The partial isometry U satisfies U^* U = P, the support projection of |A|, and is unitary on the subspace \overline{\operatorname{ran}(|A|)}.^[15] The existence of this decomposition follows from the C*-algebra structure of B(H), where A^* A is positive self-adjoint, and the spectral theorem guarantees a unique positive square root |A|.^[15] The partial isometry U is then constructed by extending an isometry from \overline{\operatorname{ran}(|A|)} to \overline{\operatorname{ran}(A)} and setting it to zero on the orthogonal complement.^[16] This ensures U preserves norms on its initial space and aligns the ranges appropriately.^[15] Uniqueness holds for |A| due to the spectral theorem, and for U on \operatorname{ran}(A), with U^* U equaling the support projection of |A|.^[15]^[16] Key properties include the equality of operator norms \|A\| = \|U\| = \||A\|\|, since U is an isometry on the relevant subspace.^[15] Spectrum relations from the matrix case extend, with the spectrum of A related to that of |A| via the partial isometry action, adjusted for possible kernel and cokernel effects in infinite dimensions.^[15] A representative example is the unilateral right shift operator S on \ell^2(\mathbb{N}), defined by S(e_n) = e_{n+1} for the standard basis \{e_n\}, which is a bounded isometry with S^* S = I but S S^* \neq I.^[15] Here, |S| = I, so the polar decomposition is S = S \cdot I, where S itself serves as the partial isometry, with initial projection I and final projection I - |e_1\rangle\langle e_1|.^[15]

Unbounded Operators

In the context of unbounded operators on a Hilbert space \mathcal{H}, the polar decomposition applies to closed densely defined linear operators A: \mathcal{D}(A) \subseteq \mathcal{H} \to \mathcal{H}, where \mathcal{D}(A) is dense in \mathcal{H}. Such an operator admits a unique factorization A = U |A|, with |A| = (A^* A)^{1/2} denoting the positive self-adjoint square root of the closure of A^* A, defined on the domain \mathcal{D}(|A|) = \mathcal{D}(A). Here, U is a partial isometry (bounded linear operator) satisfying \ker U = \ker A, with initial space \overline{\mathrm{ran}(|A|)} and final space \overline{\mathrm{ran}(A)}, such that U^* U is the orthogonal projection onto \overline{\mathrm{ran}(|A|)} and U U^* is the orthogonal projection onto \overline{\mathrm{ran}(A)}.^[17]^[18] The construction relies on the spectral theorem for the self-adjoint operator |A|, which ensures the existence of a unique positive self-adjoint extension. Specifically, since A is closed and densely defined, A^* exists and is closed, and the symmetric operator A^* A (initially defined on \mathcal{D}(A^* A) = \mathcal{D}(A) \cap \mathcal{D}(A^*)) has a self-adjoint closure, allowing |A| to be well-defined via functional calculus. The partial isometry U is then constructed by extending a densely defined operator from \mathrm{ran}(|A|^{1/2}) to the whole space: for x \in \mathcal{D}(A), set U |A|^{1/2} y = A (|A|^{-1/2} y) where y \in \mathcal{D}(|A|^{1/2}) with \mathcal{D}(A) \subseteq \mathcal{D}(|A|^{1/2}), and this extends uniquely to a bounded partial isometry on \mathcal{H}. The equality A = U |A| holds precisely on \mathcal{D}(A).^[17]^[18]^[19] Existence follows directly from the closedness of A, which implies that the closure of A^* A is self-adjoint and positive, guaranteeing the spectral resolution for |A|; the partial isometry U then extends uniquely due to the closed graph theorem applied to the graph of A. However, significant challenges arise with domains: while \mathcal{D}(|A|) = \mathcal{D}(A), the domain of U is the entire \mathcal{H}, so the product U |A| is only equal to A on \mathcal{D}(A) and may differ elsewhere, reflecting the unbounded growth of |A| outside this domain. In quantum mechanics, this manifests in operators like the momentum operator P = -i \frac{d}{dx} on L^2(\mathbb{R}) with domain the Sobolev space H^1(\mathbb{R}); its polar decomposition is P = U |P|, where in the Fourier basis, |P| multiplies by |k| and U by \mathrm{sign}(k), both preserving the domain H^1(\mathbb{R}), but illustrating how spectral measures handle unbounded spectra while respecting domain restrictions for physical observables.^[17]^[18]^[19]

Quaternion Polar Decomposition

The quaternion algebra \mathbb{H} is a four-dimensional associative division algebra over the real numbers, spanned by the basis \{1, i, j, k\} with multiplication rules i^2 = j^2 = k^2 = -1, ij = k, jk = i, and ki = j. Matrices over \mathbb{H}, denoted M_n(\mathbb{H}), are n \times n arrays with entries in \mathbb{H}, equipped with the standard Hermitian inner product \langle x, y \rangle = y^* x, where ^* denotes the conjugate transpose and the quaternion conjugate of q = a + bi + cj + dk is \overline{q} = a - bi - cj - dk. The polar decomposition for a matrix A \in M_n(\mathbb{H}) expresses A = UP, where U \in M_n(\mathbb{H}) is unitary (U^* U = I) and P \in M_n(\mathbb{H}) is positive semidefinite Hermitian (P = P^* and x^* P x \geq 0 for all x \in \mathbb{H}^n, with equality to zero only if x = 0 in the invertible case). This decomposition exists for all A \in M_n(\mathbb{H}), including non-invertible matrices, and generalizes the complex case while accounting for the non-commutativity of \mathbb{H}. The positive semidefinite factor is constructed as P = \sqrt{A^* A}, where the square root is defined via the spectral theorem for Hermitian quaternion matrices, which admit a diagonalization A^* A = V D V^* with V unitary and D real diagonal with non-negative entries, yielding \sqrt{A^* A} = V \sqrt{D} V^* using the non-negative real square root on the spectrum. The unitary factor is then U = A P^\dagger for the Moore-Penrose pseudoinverse P^\dagger in the general case, or U = A P^{-1} when A is invertible; this ensures U is a partial isometry with initial space matching the range of P. Uniqueness holds for the positive semidefinite factor P in the minimal polar decomposition (where the kernel of U is minimized), analogous to the complex setting, but non-commutativity introduces distinctions between right polar decompositions (A = UP) and left variants (A = PU'), with the standard form being right-sided. For invertible A, both U and P are unique. A simple example is the $1 \times 1 case, where any nonzero quaternion q \in \mathbb{H} decomposes as q = |q| \cdot u with |q| = \sqrt{q \overline{q}} \geq 0 the real modulus (positive semidefinite) and u = q / |q| a unit quaternion (\overline{u} u = 1), representing a scaling by a positive real followed by a rotation in 3D space.

Alternative Planar Decompositions

In the planar case, where A \in \mathbb{R}^{2 \times 2} is an invertible matrix, the polar decomposition A = UP separates A into an orthogonal factor U \in O(2) (with \det(U) = \pm 1) and a symmetric positive definite factor P, capturing rotation or reflection followed by a stretch along principal axes. Alternative factorizations exist that attempt similar separations into a rotational component and a scaling or distortion component, but they differ in structure and geometric interpretation. For instance, the QR decomposition expresses A = QR, where Q \in O(2) is orthogonal and R is upper triangular with positive diagonal entries, effectively combining anisotropic scaling along coordinate axes with shear. This contrasts with polar decomposition, as R is generally not symmetric, mixing shear and scaling effects rather than isolating pure stretch. A key distinction arises when reflections are involved: the polar decomposition accommodates \det(U) = -1 naturally, preserving the symmetric positive structure of P, whereas alternatives restricted to proper rotations in SO(2) (with \det = 1) would require adjusting the second factor, often resulting in a non-symmetric or non-positive component that loses the clean separation of orthogonal and stretch parts. The polar decomposition's uniqueness stems from its reliance on the symmetric positive square root of A^\top A, ensuring P aligns stretches with orthogonal principal directions invariant to coordinate rotations, a property not shared by basis-dependent alternatives like QR. For example, consider the simple shear matrix A = \begin{pmatrix} 1 & 1 \\ 0 & 1 \end{pmatrix}. Its polar decomposition yields P = \sqrt{A^\top A} = \begin{pmatrix} \frac{2}{\sqrt{5}} & \frac{1}{\sqrt{5}} \\ \frac{1}{\sqrt{5}} & \frac{3}{\sqrt{5}} \end{pmatrix} (symmetric with eigenvalues \frac{1 + \sqrt{5}}{2} and \frac{\sqrt{5} - 1}{2}) and U = A P^{-1}, a rotation by about 26.565 degrees; in contrast, the QR form is A = I \cdot A, where the identity is orthogonal but the upper triangular A embeds both shear and uniform scaling without distinguishing them. This highlights polar's preference for geometric applications, as its factors remain invariant under orthogonal transformations of the basis, providing a coordinate-free interpretation of deformation into rigid motion and strain.

Connections to Other Factorizations

The polar decomposition of a square invertible matrix A ∈ ℝ^{n×n} is closely related to the QR decomposition A = Q R, where Q is orthogonal and R is upper triangular with positive diagonal entries. The positive factor P in the polar decomposition A = U P is the unique symmetric positive definite square root of A^T A, while R is the upper Cholesky factor of A^T A (i.e., A^T A = R^T R). Thus, P and R are linked through the symmetric structure of A^T A, but P is symmetric whereas R is triangular; the unitary factor U = A P^{-1} differs from Q = A R^{-1} by the transformation R P^{-1}, which is generally not the identity unless A is orthogonal. This connection allows QR decompositions to be used in iterative algorithms for computing the polar decomposition, such as Newton's method, where each step orthogonalizes an approximate unitary factor via QR to ensure stability. For normal matrices, the polar decomposition aligns directly with the eigendecomposition. A normal matrix A satisfies A A^* = A^* A and admits an eigendecomposition A = V D V^, where V is unitary and D is diagonal with complex entries. The polar factors are then P = V |D| V^ (positive semidefinite, with |D| the diagonal of absolute values) and U = V (D / |D|) V^* (unitary, capturing the phases or signs in D / |D|). This separation isolates the magnitude (|λ_i|) in P and the directional component (e^{i arg(λ_i)}) in U, providing a canonical form that extends the eigendecomposition to non-normal cases via the more general polar structure; for non-normal matrices, the Jordan canonical form complicates direct alignment, as it involves non-diagonalizable blocks not separable into unitary and positive parts in the same way. Compared to the Schur decomposition A = Q T Q^*, where Q is unitary and T is upper triangular (revealing eigenvalues on the diagonal of T), the polar decomposition offers advantages in preserving symmetry for the positive component. The Schur "stretch" T is triangular and generally non-symmetric, which can obscure applications requiring symmetric positive operators (e.g., covariance modeling or metric definitions); in contrast, the polar P is explicitly symmetric positive definite, facilitating better structure exploitation in optimization and approximation tasks without additional symmetrization steps.

Numerical Computation

Algorithms for Matrices

One standard direct method for computing the polar decomposition of a matrix relies on its singular value decomposition (SVD). For a complex matrix A \in \mathbb{C}^{m \times n} with m \geq n, compute the thin SVD A = U \Sigma V^*, where U \in \mathbb{C}^{m \times n} has orthonormal columns, \Sigma \in \mathbb{R}^{n \times n} is diagonal with non-negative singular values, and V \in \mathbb{C}^{n \times n} is unitary. The right polar decomposition is then given by A = UP, where the positive semidefinite factor is P = V \Sigma V^* and the partial isometry (unitary if full rank) is U = U V^*. For the general rectangular case or m < n, the construction uses the appropriate thin SVD and pseudoinverse of P to ensure U = A P^\dagger.^[12]^[2] For the 2×2 case, closed-form expressions exist for nonsingular real matrices, leveraging traces and determinants to avoid full SVD computation. The positive semidefinite factor P = \sqrt{A^T A} can be obtained explicitly: let \tau = \operatorname{tr}(A^T A) and \delta = \det(A^T A); then P = \frac{A^T A + \sqrt{\delta} \, I}{\sqrt{\tau + 2 \sqrt{\delta}}}, assuming \tau^2 \neq 4\delta and positive definiteness. The orthogonal factor follows as U = A P^{-1}. These formulas arise from the explicit square root computation for 2×2 positive definite matrices via the Cayley-Hamilton theorem.^[20] When requiring the orthogonal factor U to be a proper rotation (i.e., \det(U) = 1 for real matrices with \det(A) > 0), adjust the signs of columns in the SVD singular vectors: specifically, if \det(U V^T) = -1, negate one column of V and the corresponding column of U to flip the determinant without altering the product U \Sigma V^*. This ensures compatibility with the special orthogonal group SO(n).^[12] The computational complexity of the SVD-based approach is O(n^3) for an n \times n square matrix using standard algorithms like bidiagonalization followed by QR iterations, making it suitable for moderate-sized problems. This method is backward stable, preserving relative accuracy even for highly ill-conditioned matrices where condition numbers exceed machine precision.^[12] For singular matrices, the decomposition extends naturally via SVD: P has rank equal to that of A, with zero singular values set to zero, and U acts as a partial isometry mapping the range of P unitarily onto the range of A. Compute U using the pseudoinverse P^\dagger restricted to the range of P, ensuring U = A P^\dagger on that subspace.^[2]

Iterative Methods and Convergence

Iterative methods for computing the polar decomposition of a matrix A typically approximate the unitary factor U and positive semidefinite factor P through successive refinements, avoiding direct computation of square roots or singular value decompositions, which can be expensive for large or ill-conditioned matrices. These methods are particularly useful for large-scale problems where direct methods are impractical, and they often exhibit quadratic convergence near the solution under suitable scaling. One prominent approach is the Newton method applied to the matrix sign function, which targets the unitary polar factor U = \sign(A), defined via the polar decomposition A = U P where P = (A^* A)^{1/2}.^[7] The unscaled Newton iteration for the sign function is given by

X_{k+1} = \frac{1}{2} (X_k + X_k^{-1}), \quad X_0 = A,

which converges to U for full-rank A, with quadratic convergence in a neighborhood of the solution provided the spectral radius condition is satisfied. To ensure global convergence and numerical stability, scaling is introduced, yielding the scaled variant

X_{k+1} = \frac{1}{2} (\mu_k X_k + \mu_k^{-1} X_k^{-*}),

where \mu_k is chosen based on norms of X_k and its inverse, such as \mu_k = (\|X_k\|_2 \|X_k^{-1}\|_2)^{1/4}; this requires at most n iterations for an n \times n matrix with optimal scaling and exhibits backward stability when inverses are computed accurately via bidiagonalization or QR factorization. Convergence is quadratic locally for full-rank matrices and linear globally, with the radius depending on \|A\|_2 < \sqrt{3} for the unscaled case to avoid divergence. A fixed-point iteration variant focuses on the positive polar factor P, derived from Newton's method for solving P^2 = A^* A. The iteration is

P_{k+1} = \frac{1}{2} (P_k + (A^* A) P_k^{-1}),

with P_0 an initial positive approximation (e.g., scaled identity); this converges linearly to P with rate depending on the condition number \kappa_2(A), and quadratic near P if started sufficiently close. Scaling by a factor \alpha > \|A^* A\|_2 / 2 ensures convergence within the unit disk, but the method can suffer instability without reformulation, motivating coupled iterations. Once P is obtained, U = A P^{-1} follows directly.^[21] The Denman-Beavers iteration provides a practical, stable alternative for simultaneous computation of U and P, reformulating the Newton square root iteration as a coupled pair to avoid explicit inverses in unstable forms:

Y_{k+1} = \frac{1}{2} (Y_k + Z_k^{-1}), \quad Z_{k+1} = \frac{1}{2} (Z_k + Y_k^{-1}),

initialized with Y_0 = [A^* A](/page/A^* A), Z_0 = I, converging to P = Y_\infty and P^{-1} = Z_\infty; then U = A Z_\infty. This achieves quadratic convergence for \|Y_0\|_2 < 1 after scaling and is numerically stable, preserving positive definiteness, particularly for normal matrices where the factors align closely with the SVD. Error analysis shows backward stability with relative error bounded by O(\epsilon \kappa_2(A)), where \epsilon is machine precision, avoiding direct square root computations that can amplify rounding errors.^[21] Recent advancements include block-coordinate methods, such as Riemannian block coordinate descent, which optimize over blocks of the unitary factor on the Stiefel manifold to compute sparse or structured polar decompositions efficiently for large-scale data; these exhibit linear convergence under restricted isometry assumptions and are suitable for high-dimensional applications like tensor approximation. Overall, these iterations offer backward stable approximations with convergence rates scaling favorably for well-conditioned problems, though ill-conditioning requires careful scaling to maintain quadratic behavior.^[22]

References

[1]
[PDF] Polar decomposition - MIT Mathematics
Nov 25, 2019 · Here is the polar decomposition I stated in class. Theorem 1.3. In the setting of (1.1), there is a factorization. T = SP. (S ∈ L(V ...
[2]
[PDF] The polar decomposition— properties, applications and algorithms
The polar decomposition was introduced by Autonne [1] in 1902. A thor- ougłi cliscussion of the history of it is given in Horn and Johnson [29, Sect. 3.0] ...Missing: theorem origin
[3]
[PDF] Lecture 7.4: Polar decomposition
Theorem. Every linear map A: X → U can be written as A = UP where P ≥ 0 and U is unitary. This is called the (left) polar decomposition of A.
[4]
[PDF] The Matrix Sign Decomposition and Its Relation to the Polar ...
Most results for the matrix sign decomposition have a counterpart for the polar decomposition A = UH, and vice versa. To illustrate this, we derive best ...<|control11|><|separator|>
[5]
https://arxiv.org/abs/2006.00841
[6]
[PDF] Polar Decomposition and the Closest Rotation - Julian Panetta
Feb 23, 2015 · A square n × n complex matrix A has two polar decompositions: the standard right polar decomposition. A = UP and the left polar decomposition ...
[7]
[PDF] Introduction to Continuum Mechanics
Next, a geometric interpretation is obtained for the polar decomposition decomposition, starting with the right polar decomposition (3.86)1. To this end ...
[8]
[PDF] Polar Decomposition of a Matrix
Apr 28, 2014 · We know from Theorem OD in a First Course in Linear Algebra that our matrices P and H are both orthonormally diagonalizable.
[9]
[PDF] 2.4.1 Stretch Vector and Stretch - ME338A CONTINUUM MECHANICS
Apr 22, 2008 · The polar decomposition theorem states that any non-singular, ... This makes the geometric interpretation of the polar decom- position ...
[10]
[PDF] 3·Singular Value Decomposition
Apr 1, 2017 · The singular value decomposition also gives an immediate formula for the 2-norm of a matrix. kAxk kxk = s1.
[11]
[PDF] Chapter 12 Singular Value Decomposition and Polar Form
Given an SVD decomposition A = VDU>, let. R = V U> and S = UDU>. It is clear that R is orthogonal and that S is positive semidefinite symmetric, and. RS = V U> ...
[12]
https://www.cis.upenn.edu/~cis5150/cis515-12-sl12.pdf
[13]
[PDF] Functional Analysis: Spectral Theory
Feb 5, 2012 · Functional Analysis: Spectral Theory. V.S. Sunder. Institute of Mathematical Sciences. Madras 600113. INDIA. July 31, 2000. Page 2. i ha ha ...Missing: textbook | Show results with:textbook
[14]
[PDF] Functional Analysis Princeton University MAT520 Lecture Notes
Aug 18, 2023 · 10 The spectral theorem for bounded normal operators. 110. 10.1 ... Functional analysis is the branch of mathematics that is obtained ...
[15]
[PDF] Chapter 1 Hilbert space and linear operators
Obviously, the following definitions and results are also valid for bounded linear operators. ... Theorem 1.6.5 (Polar decomposition). For any B ∈ B(H) there ...
[16]
Perturbation Theory for Linear Operators - SpringerLink
Free delivery 14-day returnsPerturbation Theory for Linear Operators · Reviews. "The monograph by T. Kato is an excellent textbook in the theory of linear operators in Banach and Hilbert ...
[17]
Linear Operators in Hilbert Spaces - SpringerLink
The purpose of this book is to give an introduction to the theory of linear operators on Hilbert spaces and then to proceed to the interesting applica tions of ...
[18]
https://link.springer.com/book/10.1007/978-1-4612-6027-1
[19]
[PDF] Matrix Animation and Polar Decomposition - Graphics Interface
This paper presents a method for decomposing matrices, focusing on rotation extraction using Polar Decomposition, which is useful for animation and ...
[20]
Illustration of polar decomposition | University of Utah CSM Group
Feb 14, 2013 · This posting explains the meaning of a polar decomposition, and it gives two numerical methods for computing it.
[21]
Computing the Polar Decomposition—with Applications - SIAM.org
Computing the Singular Value Decomposition of a Product of Two Matrices ... A quadratically convergent Newton method for computing the polar decomposition of a ...
[22]
[PDF] Stable and Efficient Computation of Generalized Polar ... - arXiv
Apr 14, 2021 · Abstract: We present methods for computing the generalized polar decomposition of a matrix based on the dynamically weighted Halley (DWH) ...
[23]
[PDF] arXiv:2011.12449v1 [math.NA] 25 Nov 2020
Nov 25, 2020 · ric eigendecomposition, singular value decomposition, polar decomposition, and CS ... Let B ∈ Cm×m be a nonsingular normal matrix. Let ...
[24]
The Canonical Generalized Polar Decomposition - SIAM.org
We present methods for computing the generalized polar decomposition of a matrix based on the dynamically weighted Halley iteration. This method is well ...
[25]
https://www.sciencedirect.com/science/article/pii/0024379581900239
[26]
[PDF] A Riemannian Block Coordinate Descent Method for Computing the ...
Sep 27, 2021 · the the polar decomposition as the retraction, L1 =1+. √. 2/2,L2 = √. 10. 3 A Riemannian Block Coordinate Descent Algorithm for Com- puting ...