Add KAN decomposition function

docs/src/spacetime_algebra.md

@@ -27,13 +27,13 @@ this more general algebra.
 
 One crucial new feature is the *spinor norm*: for complex quaternions,
 the GA reverse gives ``𝐐𝐐̃ = w² + x² + y² + z²``, which is a
-*"complex"* scalar, not the Euclidean norm ``|w|² + |x|² + |y|² +
-|z|²`` which will always be real.  The surprising point is that the
-normalization condition ``𝐐𝐐̃ = 1`` is now *two real* conditions,
-because this says that the "imaginary" part is zero.  This reduces the
-eight real degrees of freedom in a complex quaternion down to six,
-which is the correct number for a Lorentz transformation in four
-dimensions.
+*complex* scalar, distinct from the Euclidean norm ``|w|² + |x|² +
+|y|² + |z|²`` which will always be real.  The surprising point is that
+the normalization condition ``𝐐𝐐̃ = 1`` is now *two real*
+conditions, because this says that the "imaginary" part is zero.  This
+reduces the eight real degrees of freedom in a complex quaternion down
+to six, which is the correct number for a Lorentz transformation in
+four dimensions.
 
 
 ## The spacetime algebra
@@ -197,8 +197,6 @@ which is not equal to ``1`` in general.  Phase factors therefore do *not*
 belong to the Lorentz rotor group, even though they are unit elements under
 the Euclidean norm (``|\exp(i\varphi)| = 1``).
 
-
-
 !!! note "Real vs. complex quaternions"
 
     For real quaternions, ``𝐐\,\widetilde{𝐐} = \sum r_i^2 = \sum |r_i|^2``,
@@ -247,11 +245,177 @@ Lorentz transformations.
     and confirming that a positive rapidity boosts ``t \to \cosh\varphi\,t +
     \sinh\varphi\,x`` and ``x \to \sinh\varphi\,t + \cosh\varphi\,x``.
 
+
+## [Iwasawa's ``KA\,N`` decomposition](@id iwasawa-kan)
+
+The [Iwasawa
+decomposition](https://en.wikipedia.org/wiki/Iwasawa_decomposition) is
+a factorization of the Lorentz group or its double cover
+``\mathrm{Spin}^+(3,1)`` (or any semisimple Lie group) into three
+subgroups:
+
+```math
+G = {KA\,N}.
+```
+
+This decomposition is most useful when we have a preferred time axis
+and preferred spatial direction.  We use those to pick out boosts
+along that spatial direction, and null rotations that fix the
+corresponding null vector.  Specifically, we have the following
+subgroups:
+
+- **``K``** is the maximal compact subgroup, which in the case of the
+  Lorentz group is isomorphic to the rotation group
+  ``\mathrm{SO}(3)``, and in the case of ``\mathrm{Spin}(3,1)`` is
+  isomorphic to ``\mathrm{Spin}(3)``.
+- **``A``** is the abelian subgroup of boosts in a fixed direction.
+  We take this direction to be the z-axis, so that we have ``A =
+  \left\{ \exp\left[\tfrac{φₐ}{2} \, 𝐭𝐳 \right] \mid φₐ ∈ ℝ
+  \right\}``.
+- **``N``** is the nilpotent subgroup of null rotations, which
+  consists of all Lorentz transformations that can be represented as
+  null rotations about a fixed null vector.  We take this vector to be
+  the null vector ``\boldsymbol{ℓ}``, so that we have ``N = \left\{
+  \exp\left[\tfrac{1}{2} \boldsymbol{ℓ ξ}\right] \mid \boldsymbol{ξ} =
+  ξˣ𝐱 + ξʸ𝐲 \right\}``, where ``\boldsymbol{ℓ} = (𝐭+𝐳)/√2``.
+
+Note that ``φₐ`` here is *not* the rapidity of the overall boost if we
+factor a transformation just as a boost and a rotation.  Rather, it is
+the rapidity of the particular boost picked out as ``A`` by this
+construction.  ``N`` also contributes to the overall boost, so the
+overall rapidity is not just ``φₐ``.
+
+One nice feature of this decomposition is that we can compute the
+factorization of a given Lorentz transformation ``Λ`` fairly simply.
+We first introduce the idempotent
+
+```math
+u₊ = \frac{1}{2} (1 + 𝐭𝐳).
+```
+
+It is easy to verify ``u₊𝐭𝐳 = 𝐭𝐳u₊ = u₊`` and hence the idempotent
+identity: ``u₊² = u₊``.  Next, we write
+
+```math
+Λ = Rₖ Rₐ Rₙ,
+```
+
+where ``Rₖ ∈ K``, ``Rₐ ∈ A``, and ``Rₙ ∈ N``.  Starting with the
+rightmost factor, we have
+
+```math
+Rₙ = \exp\left[\tfrac{1}{2} \boldsymbol{ℓ ξ}\right] = 1 + \frac{1}{2} \boldsymbol{ℓ ξ},
+```
+
+because ``\boldsymbol{ℓ}`` and ``\boldsymbol{ξ}`` anticommute, and
+because ``\boldsymbol{ℓ}² = 0``, their product ``\boldsymbol{ℓ ξ}``
+itself is also nilpotent.  (Remember, that's what the "N" stands for.)
+We can simply compute the product ``Rₙ u₊`` by using the fact that
+``\boldsymbol{ξ} u₊ =  u₊ \boldsymbol{ξ}`` and expanding terms in our
+basis, then find
+
+```math
+Rₙ u₊ %&= \left(1 + \frac{1}{2} \boldsymbol{ℓ ξ}\right) u₊, \\
+%  &= u₊ + \frac{1}{2} \boldsymbol{ℓ ξ} u₊, \\
+%  &= u₊ + \frac{1}{2} \boldsymbol{ℓ} u₊ \boldsymbol{ξ}, \\
+%  &= u₊ + \frac{1}{4\sqrt{2}} (𝐭+𝐳) (1+𝐭𝐳) \boldsymbol{ξ}, \\
+%  &= u₊ + \frac{1}{4\sqrt{2}} (𝐭+𝐳+𝐭𝐭𝐳+𝐳𝐭𝐳) \boldsymbol{ξ}, \\
+%  &= u₊ + \frac{1}{4\sqrt{2}} (𝐭+𝐳-𝐳-𝐭) \boldsymbol{ξ}, \\
+= u₊.
+```
+
+Now, with
+
+```math
+Rₐ = \exp\left[\tfrac{φₐ}{2} \, 𝐭𝐳 \right]
+= \cosh\left(\frac{φₐ}{2}\right) + \sinh\left(\frac{φₐ}{2}\right) 𝐭𝐳,
+```
+
+and the fact that ``𝐭𝐳u₊ = u₊``, we have
+
+```math
+\begin{aligned}
+Rₐ u₊
+&= \left[\cosh\left(\frac{φₐ}{2}\right) + \sinh\left(\frac{φₐ}{2}\right)\right] u₊
+&= e^{φₐ/2} u₊.
+\end{aligned}
+```
+
+Putting these together, we have
+
+```math
+Λ u₊ = e^{φₐ/2} Rₖ u₊.
+```
+
+Recall that ``e^{φₐ/2}`` is a strictly positive real number, and
+``Rₖ`` is a pure rotation so it is "ℂ-real" — meaning that it is a
+linear combination of only basis elements that do not have a factor of
+``𝐭`` — whereas ``u₊`` is just a combination of 1 and a "ℂ-imaginary"
+part.  Therefore, we can take the real part of this expression to find
+
+```math
+ℂ\Re\{Λ u₊\} = \frac{e^{φₐ/2}}{2} Rₖ.
+```
+
+This is just a strictly positive real number times ``Rₖ``, which we
+know must have unit magnitude, so we can find ``Rₖ`` by normalizing
+this real part:
+
+```math
+Rₖ = \mathrm{normalize}\left(ℂ\Re\{Λ u₊\}\right).
+```
+
+We can continue solving for the remaining factors ``Rₐ`` and ``Rₙ``.
+One option would be to take the norm of the expression above and solve
+for ``φₐ``, but it requires no transcendental functions and is more
+numerically stable to extract the components directly from the
+remaining product.  We can evaluate ``RₐRₙ = Rₖ^{-1} Λ``.  But now,
+let us examine the left-hand side analytically:
+
+```math
+Rₐ Rₙ = \left[\cosh\left(\frac{φₐ}{2}\right) + \sinh\left(\frac{φₐ}{2}\right) 𝐭𝐳\right]
+  \left[1 + \frac{1}{2} \boldsymbol{ℓ ξ}\right]
+```
+
+The ``\boldsymbol{ℓ ξ}`` term is entirely within the span of ``\{𝐭𝐱,
+𝐭𝐲, 𝐳𝐱, 𝐳𝐲\}``; therefore, it is not hard to show that the
+scalar and ``𝐭𝐳`` terms in this product equal ``Rₐ``.  As usual, we
+write ``⟨·⟩₀`` to extract the scalar (grade‑0) part, and we have
+
+```math
+Rₐ = \left⟨ Rₖ^{-1} Λ \right⟩₀ + \left⟨ Rₖ^{-1} Λ\, 𝐭𝐳 \right⟩₀\, 𝐭𝐳.
+```
+
+Finally, we can find the remaining factor by rearranging the original
+factorization:
+
+```math
+Rₙ = Rₐ^{-1} Rₖ^{-1} Λ.
+```
+
+This decomposition is implemented in the [`Quaternionic.KAN`](@ref)
+function.
+
 ## API reference
 
-```@autodocs
-Modules = [Quaternionic]
-Pages   = ["Lorentz.jl"]
+```@meta
+CurrentModule = Quaternionic
+```
+
+```@docs
+Lorentz
+Lorentz(::AbstractVector)
+Boost
+ga_components
+RB
+BR
+Rv
+vR
+KAN
+ℂreal
+ℂimag
+ℂreim
+ℂconj
 ```
 
 ## Further reading

src/Lorentz.jl

@@ -268,7 +268,7 @@ This is distinct from the GA "real" part, given by `real(Λ)`, which is construc
 reverse, rather than the complex-conjugate.  That will just be the scalar part of `Λ`, but
 may be complex; this function returns a full quaternion, with all components real.
 """
-ℂreal(Λ::Lorentz{T}) where {T<:Real} = Quaternion{T}(
+ℂreal(Λ::AbstractQuaternion{Complex{T}}) where {T<:Real} = Quaternion{T}(
     real(Λ[1]), real(Λ[2]), real(Λ[3]), real(Λ[4])
 )
 
@@ -285,7 +285,7 @@ the reverse, rather than the complex-conjugate.  That will just be the quaternio
 part of `Λ`, but may be complex; this function returns a full quaternion, with all
 components real.
 """
-ℂimag(Λ::Lorentz{T}) where {T<:Real} = Quaternion{T}(
+ℂimag(Λ::AbstractQuaternion{Complex{T}}) where {T<:Real} = Quaternion{T}(
     imag(Λ[1]), imag(Λ[2]), imag(Λ[3]), imag(Λ[4])
 )
 
@@ -323,7 +323,8 @@ find `cosh(η/2)` as the magnitude of the complex-real part of `Λ`.  Then, we f
 `sinh(η/2)*n̂` simply by multiplying the complex-imaginary part of `Λ` by `inv(R)`, and
 reconstruct `B` simply by adding the result to `cosh(η/2)`.
 
-See also [`BR`](@ref).
+See also [`BR`](@ref), [`vR`](@ref), and [`Rv`](@ref) for similar forms, and [`KAN`](@ref)
+for the Iwasawa decomposition.
 """
 function RB(Λ::Lorentz{T}) where {T<:Real}
     ℜΛ, ℑΛ = ℂreim(Λ)
@@ -339,6 +340,9 @@ Polar decomposition `Λ = B * R`: pure boost `B` followed by pure rotation `R`.
 
 The algorithm here is the same as for [`RB`](@ref), but with the order of multiplication by
 `inv(R)` reversed.
+
+See also [`vR`](@ref) and [`Rv`](@ref) for similar forms, and [`KAN`](@ref) for the Iwasawa
+decomposition.
 """
 function BR(Λ::Lorentz{T}) where {T<:Real}
     ℜΛ, ℑΛ = ℂreim(Λ)
@@ -370,7 +374,8 @@ where the last equality uses the double-angle identity ``\cosh(η) = 2\cosh^2(η
 The last form is made up of entirely of the scalar component of `B`, does not involve any
 cancellation or division by small numbers, and avoids computing the norm of the vector part.
 
-See also [`vR`](@ref), [`RB`](@ref), and [`BR`](@ref).
+See also [`vR`](@ref), [`RB`](@ref), and [`BR`](@ref) for similar forms, and [`KAN`](@ref)
+for the Iwasawa decomposition.
 """
 function Rv(Λ::Lorentz{T}) where {T<:Real}
     R, B = RB(Λ)
@@ -386,7 +391,8 @@ end
 Return the boost velocity `v⃗` and pure rotation `R` such that `Λ = Boost(η, v̂) * R`.  See
 [`Rv`](@ref) for more details.
 
-See also [`BR`](@ref) and [`RB`](@ref).
+See also [`BR`](@ref) and [`RB`](@ref) for similar forms, and [`KAN`](@ref) for the Iwasawa
+decomposition.
 """
 function vR(Λ::Lorentz{T}) where {T<:Real}
     B, R = BR(Λ)
@@ -395,3 +401,45 @@ function vR(Λ::Lorentz{T}) where {T<:Real}
     v⃗ = v̂sinhη╱2 * (2coshη╱2 / (2coshη╱2^2 - 1))
     return v⃗, R
 end
+
+# ---------------------------------------------------------------------------
+# Iwasawa's KAN decomposition
+# ---------------------------------------------------------------------------
+
+@doc raw"""
+    KAN(Λ::Lorentz{T}) → (Rₖ, Rₐ, Rₙ)
+
+Compute the Iwasawa ``KAN`` decomposition of the Lorentz transformation `Λ`, returning the
+three factors such that `Λ = Rₖ * Rₐ * Rₙ`, where
+
+- `Rₖ ∈ K` is a pure rotation (the maximal compact factor),
+- `Rₐ ∈ A` is a boost along the ``𝐳`` axis, `Rₐ = cosh(φₐ/2) + sinh(φₐ/2) 𝐭𝐳`, and
+- `Rₙ ∈ N` is a null rotation fixing the null vector ``ℓ = (𝐭+𝐳)/√2``.
+
+See the [`KAN` decomposition section](@ref iwasawa-kan) of the spacetime-algebra
+documentation for a full discussion and derivation of the method used in this function.
+
+See also [`BR`](@ref), [`RB`](@ref), [`vR`](@ref), and [`Rv`](@ref) for polar
+decompositions.
+"""
+function KAN(Λ::Lorentz{T}) where {T<:Real}
+    𝐭𝐳 = im * 𝐤
+    u₊ = (1 + 𝐭𝐳) / T(2)
+    # We readily find the rotation factor by projecting with the idempotent, taking the pure
+    # spatial part, and normalizing.
+    ℂℜΛu₊ = ℂreal(Λ * u₊)
+    Rₖ = rotor(ℂℜΛu₊)
+    # We avoid transcendental functions `ln`, `sinh`, and `cosh` by using the structure of
+    # the following product to read off just the relevant components.
+    RₐRₙ = conj(Rₖ) * Λ
+    coshφ╱2, sinhφ╱2 = real(RₐRₙ.w), imag(RₐRₙ.z)
+    Rₐ = coshφ╱2 + sinhφ╱2 * 𝐭𝐳
+    # Rather than returning the raw `Rₐ⁻¹ * RₐRₙ`, which may have accumulated roundoff
+    # errors, we project that result onto the null-rotation sector by averaging 𝐢 and 𝐣
+    # components, and setting 𝟏 and 𝐤 components to their exact values.
+    Rₙraw = conj(Rₐ) * RₐRₙ
+    ξˣ╱2sqrt2 = (imag(Rₙraw.x) - real(Rₙraw.y)) / 2
+    ξʸ╱2sqrt2 = (imag(Rₙraw.y) + real(Rₙraw.x)) / 2
+    Rₙ = Lorentz{T}(1, ξʸ╱2sqrt2 + im * ξˣ╱2sqrt2, -ξˣ╱2sqrt2 + im * ξʸ╱2sqrt2, 0)
+    return Rₖ, Rₐ, Rₙ
+end

src/Quaternionic.jl

@@ -15,7 +15,7 @@ export QuatVec, quatvec, QuatVecF64, QuatVecF32, QuatVecF16
 export components, basetype
 public value, iszerovalue
 public ℂconj, ℂreal, ℂimag, ℂreim
-public RB, BR, Rv, vR
+public RB, BR, Rv, vR, KAN
 export (⋅), (×), (×̂), normalize, norm
 export abs2vec, absvec
 export from_float_array, to_float_array,