From 03ebd9f0e291fca71213a73283fa314dc611557d Mon Sep 17 00:00:00 2001
From: Ronny <git@ronnybergmann.net>
Date: Sat, 6 Mar 2021 17:05:22 +0100
Subject: [PATCH 01/32] reorganize code and sketch initial algorithms of exp
 and log for canonical metric on manifolds.

---
 src/Manifolds.jl                           |   1 +
 src/manifolds/Stiefel.jl                   | 140 +------------------
 src/manifolds/StiefelCanonicalMetric.jl    |  95 +++++++++++++
 src/manifolds/StiefelEuclideanMetric.jl    | 148 +++++++++++++++++++++
 src/manifolds/SymmetricPositiveDefinite.jl |   2 +-
 5 files changed, 246 insertions(+), 140 deletions(-)
 create mode 100644 src/manifolds/StiefelCanonicalMetric.jl
 create mode 100644 src/manifolds/StiefelEuclideanMetric.jl

diff --git a/src/Manifolds.jl b/src/Manifolds.jl
index 116a052a88..e29788f356 100644
--- a/src/Manifolds.jl
+++ b/src/Manifolds.jl
@@ -160,6 +160,7 @@ include("manifolds/Rotations.jl")
 include("manifolds/SkewSymmetric.jl")
 include("manifolds/Spectrahedron.jl")
 include("manifolds/Stiefel.jl")
+include("manifolds/StiefelEuclideanMetric.jl")
 include("manifolds/Sphere.jl")
 include("manifolds/SphereSymmetricMatrices.jl")
 include("manifolds/Symmetric.jl")
diff --git a/src/manifolds/Stiefel.jl b/src/manifolds/Stiefel.jl
index 243d2f2e5b..4ed2b960e2 100644
--- a/src/manifolds/Stiefel.jl
+++ b/src/manifolds/Stiefel.jl
@@ -123,105 +123,6 @@ end
 
 decorated_manifold(::Stiefel{N,K,𝔽}) where {N,K,𝔽} = Euclidean(N, K; field=𝔽)
 
-@doc raw"""
-    exp(M::Stiefel, p, X)
-
-Compute the exponential map on the [`Stiefel`](@ref)`{n,k,𝔽}`() manifold `M`
-emanating from `p` in tangent direction `X`.
-
-````math
-\exp_p X = \begin{pmatrix}
-   p\\X
- \end{pmatrix}
- \operatorname{Exp}
- \left(
- \begin{pmatrix} p^{\mathrm{H}}X & - X^{\mathrm{H}}X\\
- I_n & p^{\mathrm{H}}X\end{pmatrix}
- \right)
-\begin{pmatrix}  \exp( -p^{\mathrm{H}}X) \\ 0_n\end{pmatrix},
-````
-
-where $\operatorname{Exp}$ denotes matrix exponential,
-$\cdot^{\mathrm{H}}$ denotes the complex conjugate transpose or Hermitian, and $I_k$ and
-$0_k$ are the identity matrix and the zero matrix of dimension $k × k$, respectively.
-"""
-exp(::Stiefel, ::Any...)
-
-function exp!(::Stiefel{n,k}, q, p, X) where {n,k}
-    return copyto!(
-        q,
-        [p X] *
-        exp([p'X -X'*X; one(zeros(eltype(p), k, k)) p'*X]) *
-        [exp(-p'X); zeros(eltype(p), k, k)],
-    )
-end
-
-@doc raw"""
-    get_basis(M::Stiefel{n,k,ℝ}, p, B::DefaultOrthonormalBasis) where {n,k}
-
-Create the default basis using the parametrization for any $X ∈ T_p\mathcal M$.
-Set $p_\bot \in ℝ^{n\times(n-k)}$ the matrix such that the $n\times n$ matrix of the common
-columns $[p\ p_\bot]$ is an ONB.
-For any skew symmetric matrix $a ∈ ℝ^{k\times k}$ and any $b ∈ ℝ^{(n-k)\times k}$ the matrix
-
-````math
-X = pa + p_\bot b ∈ T_p\mathcal M
-````
-
-and we can use the $\frac{1}{2}k(k-1) + (n-k)k = nk-\frac{1}{2}k(k+1)$ entries
-of $a$ and $b$ to specify a basis for the tangent space.
-using unit vectors for constructing both
-the upper matrix of $a$ to build a skew symmetric matrix and the matrix b, the default
-basis is constructed.
-
-Since $[p\ p_\bot]$ is an automorphism on $ℝ^{n\times p}$ the elements of $a$ and $b$ are
-orthonormal coordinates for the tangent space. To be precise exactly one element in the upper
-trangular entries of $a$ is set to $1$ its symmetric entry to $-1$ and we normalize with
-the factor $\frac{1}{\sqrt{2}}$ and for $b$ one can just use unit vectors reshaped to a matrix
-to obtain orthonormal set of parameters.
-"""
-function get_basis(M::Stiefel{n,k,ℝ}, p, B::DefaultOrthonormalBasis{ℝ}) where {n,k}
-    V = get_vectors(M, p, B)
-    return CachedBasis(B, V)
-end
-
-function get_coordinates!(
-    M::Stiefel{n,k,ℝ},
-    c,
-    p,
-    X,
-    B::DefaultOrthonormalBasis{ℝ},
-) where {n,k}
-    V = get_vectors(M, p, B)
-    c .= [inner(M, p, v, X) for v in V]
-    return c
-end
-
-function get_vector!(M::Stiefel{n,k,ℝ}, X, p, c, B::DefaultOrthonormalBasis{ℝ}) where {n,k}
-    V = get_vectors(M, p, B)
-    zero_tangent_vector!(M, X, p)
-    length(c) < length(V) && error(
-        "Coordinate vector too short. Excpected $(length(V)), but only got $(length(c)) entries.",
-    )
-    @inbounds for i in 1:length(V)
-        X .+= c[i] .* V[i]
-    end
-    return X
-end
-
-function get_vectors(::Stiefel{n,k,ℝ}, p, ::DefaultOrthonormalBasis{ℝ}) where {n,k}
-    p⊥ = nullspace([p zeros(n, n - k)])
-    an = div(k * (k - 1), 2)
-    bn = (n - k) * k
-    V = vcat(
-        [p * vec2skew(1 / sqrt(2) .* _euclidean_unit_vector(an, i), k) for i in 1:an],
-        [p⊥ * reshape(_euclidean_unit_vector(bn, j), (n - k, k)) for j in 1:bn],
-    )
-    return V
-end
-
-_euclidean_unit_vector(n, i) = [k == i ? 1.0 : 0.0 for k in 1:n]
-
 @doc raw"""
     inverse_retract(M::Stiefel, p, q, ::PolarInverseRetraction)
 
@@ -381,42 +282,6 @@ manifold_dimension(::Stiefel{n,k,ℝ}) where {n,k} = n * k - div(k * (k + 1), 2)
 manifold_dimension(::Stiefel{n,k,ℂ}) where {n,k} = 2 * n * k - k * k
 manifold_dimension(::Stiefel{n,k,ℍ}) where {n,k} = 4 * n * k - k * (2k - 1)
 
-@doc raw"""
-    project(M::Stiefel,p)
-
-Projects `p` from the embedding onto the [`Stiefel`](@ref) `M`, i.e. compute `q`
-as the polar decomposition of $p$ such that $q^{\mathrm{H}q$ is the identity,
-where $\cdot^{\mathrm{H}}$ denotes the hermitian, i.e. complex conjugate transposed.
-"""
-project(::Stiefel, ::Any, ::Any)
-
-function project!(::Stiefel, q, p)
-    s = svd(p)
-    mul!(q, s.U, s.Vt)
-    return q
-end
-
-@doc raw"""
-    project(M::Stiefel, p, X)
-
-Project `X` onto the tangent space of `p` to the [`Stiefel`](@ref) manifold `M`.
-The formula reads
-
-````math
-\operatorname{proj}_{\mathcal M}(p, X) = X - p \operatorname{Sym}(p^{\mathrm{H}}X),
-````
-
-where $\operatorname{Sym}(q)$ is the symmetrization of $q$, e.g. by
-$\operatorname{Sym}(q) = \frac{q^{\mathrm{H}}+q}{2}$.
-"""
-project(::Stiefel, ::Any...)
-
-function project!(::Stiefel, Y, p, X)
-    A = p' * X
-    copyto!(Y, X - p * Hermitian((A + A') / 2))
-    return Y
-end
-
 @doc raw"""
     retract(::Stiefel, p, X, ::CayleyRetraction)
 
@@ -812,7 +677,4 @@ function vector_transport_to!(
         Y,
         q * (UpperTriangular(qtXrf) - UpperTriangular(qtXrf)') + Xrf - q * qtXrf,
     )
-end
-function vector_transport_to!(M::Stiefel, Y, ::Any, X, q, ::ProjectionTransport)
-    return project!(M, Y, q, X)
-end
+end
\ No newline at end of file
diff --git a/src/manifolds/StiefelCanonicalMetric.jl b/src/manifolds/StiefelCanonicalMetric.jl
new file mode 100644
index 0000000000..b5f2b93ada
--- /dev/null
+++ b/src/manifolds/StiefelCanonicalMetric.jl
@@ -0,0 +1,95 @@
+@doc raw"""
+    CanonicalMetric <: Metric
+
+The Canonical Metric is name used on several manifolds, for example for the [`Stiefel`](@ref)
+manifold, see for example[^EdelmanAriasSmith1998].
+
+[^EdelmanAriasSmith1998]:
+    > Edelman, A., Ariar, T. A., Smith, S. T.:
+    > _The Geometry of Algorihthms with Orthogonality Constraints_,
+    > SIAM Journal on Matrix Analysis and Applications (20(2), pp. 303–353, 1998.
+    > doi: [10.1137/S0895479895290954](https://doi.org/10.1137/S0895479895290954)
+    > arxiv: [9806030](https://arxiv.org/abs/physics/9806030)
+"""
+struct CanonicalMetric <: RiemannianMetric end
+
+"""
+    q = exp(MetricManifold{ℝ, Stiefel{n,k,ℝ}, CanonicalMetric}, p, X)
+    exp!(MetricManifold{ℝ, Stiefel{n,k,ℝ}, q, CanonicalMetric}, p, X)
+
+Compute the exponential map on the [`Stiefel`](@ref)`(n,k)` manifold with respect to the [`CanonicalMetric`](@ref).
+
+First, decompose The tangent vector ``X`` into its horizontal and vertical component with
+respect to ``p``, i.e.
+
+```math
+X = pp^{\mathrm{T}}X + (I_n-pp^{\mathrm{T}})X,
+```
+where ``I_n`` is the ``n\times n`` identity matrix.
+We introduce ``A=p^{\mathrm{T}}X`` and ``QR = (I_n-pp^{\mathrm{T}})X``` the qr decomposition
+of the vertical component. Then using the matrix exponential ``operatorname{Exp}`` we introduce ``B`` and ``C`` as
+
+```math
+\begin{pmatrix}
+B\\C
+\end{pmatrix}
+\coloneqq
+\operatorname{Exp}\left(
+\begin{pmatrix}
+A & -R^{\mathrm{T}}\\ R & 0
+\end{pmatrix}
+\right)
+\begin{pmatrix}I_k\\0\end{pmatrix}
+```
+
+the exponential map reads
+
+```math
+q = \exp_p X = pC + QB.
+```
+For more details, see [^EdelmanAriasSmith1998][^Zimmermann2017].
+"""
+exp(::MetricManifold{ℝ, Stiefel{n,k,ℝ}, CanonicalMetric}, ::Any...) where {n,k}
+
+function exp!(::MetricManifold{ℝ, Stiefel{n,k,ℝ}, CanonicalMetric}, Q, p, X) where{n,k}
+    A = p'*X
+    QR = qr(X-p*A)
+    BC_ext = exp([A -QR.R'; QR.R, zeros(k,k)])
+    q .= [p Matrix(QR.Q)] * BC_ext[:,1:k]
+    return q
+end
+
+"""
+    X = log(MetricManifold{ℝ, Stiefel{n,k,ℝ}, CanonicalMetric}, p, q)
+    log!(MetricManifold{ℝ, Stiefel{n,k,ℝ}, X, CanonicalMetric}, p, q)
+
+Compute the logarithmic map on the [`Stiefel`](@ref)`(n,k)` manifold with respect to the [`CanonicalMetric`](@ref)
+using a matrix-algebraic based approach to an iterative inversion of the formula of the
+[`exp`](@ref exp(::MetricManifold{ℝ, Stiefel{n,k,ℝ}, CanonicalMetric}, ::Any...)).
+
+The algorithm is derived in[^Zimmermann2017].
+
+[^Zimmermann2017]:
+    > Zimmermann, R.: _A matrix-algebraic algorithm for the Riemannian logarithm on the Stiefel manifold under the canoncial metric.
+    > SIAM Journal on Matrix Analysis and Applications 28(2), pp. 322-342, 2017.
+    > doi: [10.1137/16M1074485](https://doi.org/10.1137/16M1074485),
+    > arXiv: [1604.05054](https://arxiv.org/abs/1604.05054).
+"""
+log(::MetricManifold{ℝ, Stiefel{n,k,ℝ}, CanonicalMetric}, ::Any...) where {n,k}
+
+function log!(::MetricManifold{ℝ, Stiefel{n,k,ℝ}, CanonicalMetric}, p, q; tolerance=1e-9, maxiter=1e5) where{n,k}
+    M = p'*q
+    QR = qr(q-p*M)
+    V = qr([M;Matrix(QR.R)]).Q*Matrix(I, 2*k,2*k)
+    S = svd(V[(k+1):2*k, (k+1):2k]) #bottom right corner
+    V[:,(k+1):2*k] = V[:,(k+1):2*k]*(S.V*S.U')
+    LV = log(V)
+    C = view(LV, (k+1):2*k, (k+1):2*k)
+    i=0
+    while (i < maxiter) && (norm(C) > tolerance)
+        LV = log(V)
+        V[:, (k+1):2*k] *= exp(-C)
+    end
+    X .= p*LV[1:k,1:k] + QR.Q*LV[(k+1):2*k, 1:k]
+    return X
+end
diff --git a/src/manifolds/StiefelEuclideanMetric.jl b/src/manifolds/StiefelEuclideanMetric.jl
new file mode 100644
index 0000000000..c49ba01b88
--- /dev/null
+++ b/src/manifolds/StiefelEuclideanMetric.jl
@@ -0,0 +1,148 @@
+default_metric_dispatch(::Stiefel, ::EuclideanMetric) = Val(true)
+
+function decorator_transparent_dispatch(
+    ::typeof(distance),
+    ::Stiefel,
+    args...,
+)
+    return Val(:intransparent)
+end
+
+@doc raw"""
+    exp(M::Stiefel, p, X)
+
+Compute the exponential map on the [`Stiefel`](@ref)`{n,k,𝔽}`() manifold `M`
+emanating from `p` in tangent direction `X`.
+
+````math
+\exp_p X = \begin{pmatrix}
+   p\\X
+ \end{pmatrix}
+ \operatorname{Exp}
+ \left(
+ \begin{pmatrix} p^{\mathrm{H}}X & - X^{\mathrm{H}}X\\
+ I_n & p^{\mathrm{H}}X\end{pmatrix}
+ \right)
+\begin{pmatrix}  \exp( -p^{\mathrm{H}}X) \\ 0_n\end{pmatrix},
+````
+
+where $\operatorname{Exp}$ denotes matrix exponential,
+$\cdot^{\mathrm{H}}$ denotes the complex conjugate transpose or Hermitian, and $I_k$ and
+$0_k$ are the identity matrix and the zero matrix of dimension $k × k$, respectively.
+"""
+exp(::Stiefel, ::Any...)
+
+function exp!(::Stiefel{n,k}, q, p, X) where {n,k}
+    return copyto!(
+        q,
+        [p X] *
+        exp([p'X -X'*X; one(zeros(eltype(p), k, k)) p'*X]) *
+        [exp(-p'X); zeros(eltype(p), k, k)],
+    )
+end
+
+@doc raw"""
+    get_basis(M::Stiefel{n,k,ℝ}, p, B::DefaultOrthonormalBasis) where {n,k}
+
+Create the default basis using the parametrization for any $X ∈ T_p\mathcal M$.
+Set $p_\bot \in ℝ^{n\times(n-k)}$ the matrix such that the $n\times n$ matrix of the common
+columns $[p\ p_\bot]$ is an ONB.
+For any skew symmetric matrix $a ∈ ℝ^{k\times k}$ and any $b ∈ ℝ^{(n-k)\times k}$ the matrix
+
+````math
+X = pa + p_\bot b ∈ T_p\mathcal M
+````
+
+and we can use the $\frac{1}{2}k(k-1) + (n-k)k = nk-\frac{1}{2}k(k+1)$ entries
+of $a$ and $b$ to specify a basis for the tangent space.
+using unit vectors for constructing both
+the upper matrix of $a$ to build a skew symmetric matrix and the matrix b, the default
+basis is constructed.
+
+Since $[p\ p_\bot]$ is an automorphism on $ℝ^{n\times p}$ the elements of $a$ and $b$ are
+orthonormal coordinates for the tangent space. To be precise exactly one element in the upper
+trangular entries of $a$ is set to $1$ its symmetric entry to $-1$ and we normalize with
+the factor $\frac{1}{\sqrt{2}}$ and for $b$ one can just use unit vectors reshaped to a matrix
+to obtain orthonormal set of parameters.
+"""
+function get_basis(M::Stiefel{n,k,ℝ}, p, B::DefaultOrthonormalBasis{ℝ}) where {n,k}
+    V = get_vectors(M, p, B)
+    return CachedBasis(B, V)
+end
+
+function get_coordinates!(
+    M::Stiefel{n,k,ℝ},
+    c,
+    p,
+    X,
+    B::DefaultOrthonormalBasis{ℝ},
+) where {n,k}
+    V = get_vectors(M, p, B)
+    c .= [inner(M, p, v, X) for v in V]
+    return c
+end
+
+function get_vector!(M::Stiefel{n,k,ℝ}, X, p, c, B::DefaultOrthonormalBasis{ℝ}) where {n,k}
+    V = get_vectors(M, p, B)
+    zero_tangent_vector!(M, X, p)
+    length(c) < length(V) && error(
+        "Coordinate vector too short. Excpected $(length(V)), but only got $(length(c)) entries.",
+    )
+    @inbounds for i in 1:length(V)
+        X .+= c[i] .* V[i]
+    end
+    return X
+end
+
+function get_vectors(::Stiefel{n,k,ℝ}, p, ::DefaultOrthonormalBasis{ℝ}) where {n,k}
+    p⊥ = nullspace([p zeros(n, n - k)])
+    an = div(k * (k - 1), 2)
+    bn = (n - k) * k
+    V = vcat(
+        [p * vec2skew(1 / sqrt(2) .* _euclidean_unit_vector(an, i), k) for i in 1:an],
+        [p⊥ * reshape(_euclidean_unit_vector(bn, j), (n - k, k)) for j in 1:bn],
+    )
+    return V
+end
+
+_euclidean_unit_vector(n, i) = [k == i ? 1.0 : 0.0 for k in 1:n]
+
+@doc raw"""
+    project(M::Stiefel,p)
+
+Projects `p` from the embedding onto the [`Stiefel`](@ref) `M`, i.e. compute `q`
+as the polar decomposition of $p$ such that $q^{\mathrm{H}q$ is the identity,
+where $\cdot^{\mathrm{H}}$ denotes the hermitian, i.e. complex conjugate transposed.
+"""
+project(::Stiefel, ::Any, ::Any)
+
+function project!(::Stiefel, q, p)
+    s = svd(p)
+    mul!(q, s.U, s.Vt)
+    return q
+end
+
+@doc raw"""
+    project(M::Stiefel, p, X)
+
+Project `X` onto the tangent space of `p` to the [`Stiefel`](@ref) manifold `M`.
+The formula reads
+
+````math
+\operatorname{proj}_{\mathcal M}(p, X) = X - p \operatorname{Sym}(p^{\mathrm{H}}X),
+````
+
+where $\operatorname{Sym}(q)$ is the symmetrization of $q$, e.g. by
+$\operatorname{Sym}(q) = \frac{q^{\mathrm{H}}+q}{2}$.
+"""
+project(::Stiefel, ::Any...)
+
+function project!(::Stiefel, Y, p, X)
+    A = p' * X
+    copyto!(Y, X - p * Hermitian((A + A') / 2))
+    return Y
+end
+
+function vector_transport_to!(M::Stiefel, Y, ::Any, X, q, ::ProjectionTransport)
+    return project!(M, Y, q, X)
+end
diff --git a/src/manifolds/SymmetricPositiveDefinite.jl b/src/manifolds/SymmetricPositiveDefinite.jl
index 8828f75a08..515f00364d 100644
--- a/src/manifolds/SymmetricPositiveDefinite.jl
+++ b/src/manifolds/SymmetricPositiveDefinite.jl
@@ -170,4 +170,4 @@ definite matrix `x` on the [`SymmetricPositiveDefinite`](@ref) manifold `M`.
 """
 zero_tangent_vector(::SymmetricPositiveDefinite, ::Any)
 
-zero_tangent_vector!(M::SymmetricPositiveDefinite{N}, v, x) where {N} = fill!(v, 0)
+zero_tangent_vector!(::SymmetricPositiveDefinite{N}, v, x) where {N} = fill!(v, 0)

From 1908eb697136dbfb4de90cf0f68909cdfd5e5ee4 Mon Sep 17 00:00:00 2001
From: Ronny <git@ronnybergmann.net>
Date: Sat, 6 Mar 2021 18:30:14 +0100
Subject: [PATCH 02/32] update documentation.

---
 docs/make.jl                            |  2 +-
 docs/src/manifolds/stiefel.md           | 38 ++++++++++++++
 src/manifolds/Stiefel.jl                |  2 +-
 src/manifolds/StiefelCanonicalMetric.jl | 66 ++++++++++++++++---------
 src/manifolds/StiefelEuclideanMetric.jl |  6 +--
 5 files changed, 84 insertions(+), 30 deletions(-)

diff --git a/docs/make.jl b/docs/make.jl
index c4aa2186dd..0319e66897 100644
--- a/docs/make.jl
+++ b/docs/make.jl
@@ -46,7 +46,7 @@ makedocs(
                 "Rotations" => "manifolds/rotations.md",
                 "Skew-symmetric matrices" => "manifolds/skewsymmetric.md",
                 "Spectrahedron" => "manifolds/spectrahedron.md",
-                "Sphere" => "manifolds/sphere.md",
+#                "Sphere" => "manifolds/sphere.md",
                 "Stiefel" => "manifolds/stiefel.md",
                 "Symmetric matrices" => "manifolds/symmetric.md",
                 "Symmetric positive definite" =>
diff --git a/docs/src/manifolds/stiefel.md b/docs/src/manifolds/stiefel.md
index e1312f551b..e07d10d958 100644
--- a/docs/src/manifolds/stiefel.md
+++ b/docs/src/manifolds/stiefel.md
@@ -1,9 +1,47 @@
 # Stiefel
 
+## Common and metric independent functions
+
 ```@autodocs
 Modules = [Manifolds]
 Pages = ["manifolds/Stiefel.jl"]
 Order = [:type, :function]
 ```
 
+## Default metric: the Euclidean metric
+
+The [`EuclideanMetric`](@ref) is obtained from the embedding of the Stiefel manifold in ``ℝ^{n,k}``.
+
+```@autodocs
+Modules = [Manifolds]
+Pages = ["manifolds/StiefelEuclideanMetric.jl"]
+Order = [:function]
+```
+
+## The canonical metric
+
+```@autodocs
+Modules = [Manifolds]
+Pages = ["manifolds/StiefelCanonicalMetric.jl"]
+Order = [:type]
+```
+
+Any ``X∈T_p\mathcal M``, ``p∈\mathcal M``, can be written as
+
+```math
+X = pA + (I_n-pp^{\mathrm{T}})B,
+\quad
+A ∈ ℝ^{p\imes p} \text{skew-symmetric},
+B ∈ ℝ^{n\times p} \text{arbitrary.}
+```
+
+In the Euclidean case, the elements from ``A`` are counted twice (i.e. weighted with a factor of 2).
+The canonical metric avoids this.
+
+```@autodocs
+Modules = [Manifolds]
+Pages = ["manifolds/StiefelCanonicalMetric.jl"]
+Order = [:function]
+```
+
 ## Literature
diff --git a/src/manifolds/Stiefel.jl b/src/manifolds/Stiefel.jl
index 4ed2b960e2..a59e65ff8f 100644
--- a/src/manifolds/Stiefel.jl
+++ b/src/manifolds/Stiefel.jl
@@ -677,4 +677,4 @@ function vector_transport_to!(
         Y,
         q * (UpperTriangular(qtXrf) - UpperTriangular(qtXrf)') + Xrf - q * qtXrf,
     )
-end
\ No newline at end of file
+end
diff --git a/src/manifolds/StiefelCanonicalMetric.jl b/src/manifolds/StiefelCanonicalMetric.jl
index b5f2b93ada..aa6d1d1d10 100644
--- a/src/manifolds/StiefelCanonicalMetric.jl
+++ b/src/manifolds/StiefelCanonicalMetric.jl
@@ -13,9 +13,9 @@ manifold, see for example[^EdelmanAriasSmith1998].
 """
 struct CanonicalMetric <: RiemannianMetric end
 
-"""
-    q = exp(MetricManifold{ℝ, Stiefel{n,k,ℝ}, CanonicalMetric}, p, X)
-    exp!(MetricManifold{ℝ, Stiefel{n,k,ℝ}, q, CanonicalMetric}, p, X)
+@doc raw"""
+    q = exp(M::MetricManifold{ℝ, Stiefel{n,k,ℝ}, CanonicalMetric}, p, X)
+    exp!(M::MetricManifold{ℝ, Stiefel{n,k,ℝ}, q, CanonicalMetric}, p, X)
 
 Compute the exponential map on the [`Stiefel`](@ref)`(n,k)` manifold with respect to the [`CanonicalMetric`](@ref).
 
@@ -49,19 +49,33 @@ q = \exp_p X = pC + QB.
 ```
 For more details, see [^EdelmanAriasSmith1998][^Zimmermann2017].
 """
-exp(::MetricManifold{ℝ, Stiefel{n,k,ℝ}, CanonicalMetric}, ::Any...) where {n,k}
+exp(::MetricManifold{ℝ,Stiefel{n,k,ℝ},CanonicalMetric}, ::Any...) where {n,k}
 
-function exp!(::MetricManifold{ℝ, Stiefel{n,k,ℝ}, CanonicalMetric}, Q, p, X) where{n,k}
-    A = p'*X
-    QR = qr(X-p*A)
-    BC_ext = exp([A -QR.R'; QR.R, zeros(k,k)])
-    q .= [p Matrix(QR.Q)] * BC_ext[:,1:k]
+function exp!(::MetricManifold{ℝ,Stiefel{n,k,ℝ},CanonicalMetric}, Q, p, X) where {n,k}
+    A = p' * X
+    QR = qr(X - p * A)
+    BC_ext = exp([A -QR.R'; QR.R, zeros(k, k)])
+    q .= [p Matrix(QR.Q)] * BC_ext[:, 1:k]
     return q
 end
 
+@doc raw"""
+    inner(M::MetricManifold{ℝ, Stiefel{n,k,ℝ}, X, CanonicalMetric}, p,, X, Y)
+
+compute the inner procuct on the [`Stiefel`](@ref) manifold with respect to the
+[`CanonicalMetric`](@ref). The formula reads
+
+```
+g_p(X,Y) = \operatorname{tr}\bigl( X^{\mathrm{T}}(I_n - \frac{1}{2}pp^{\mathrm{T}})Y \bigr).
+```
 """
-    X = log(MetricManifold{ℝ, Stiefel{n,k,ℝ}, CanonicalMetric}, p, q)
-    log!(MetricManifold{ℝ, Stiefel{n,k,ℝ}, X, CanonicalMetric}, p, q)
+function inner(::MetricManifold{ℝ,Stiefel{n,k,ℝ},X,CanonicalMetric}, p, X, Y) where {n,k}
+    return tr(X' * (Matrix(I, n, n) - 0.5 * p * p') * Y)
+end
+
+@doc raw"""
+    X = log(M::MetricManifold{ℝ, Stiefel{n,k,ℝ}, CanonicalMetric}, p, q)
+    log!(M::MetricManifold{ℝ, Stiefel{n,k,ℝ}, X, CanonicalMetric}, p, q)
 
 Compute the logarithmic map on the [`Stiefel`](@ref)`(n,k)` manifold with respect to the [`CanonicalMetric`](@ref)
 using a matrix-algebraic based approach to an iterative inversion of the formula of the
@@ -75,21 +89,27 @@ The algorithm is derived in[^Zimmermann2017].
     > doi: [10.1137/16M1074485](https://doi.org/10.1137/16M1074485),
     > arXiv: [1604.05054](https://arxiv.org/abs/1604.05054).
 """
-log(::MetricManifold{ℝ, Stiefel{n,k,ℝ}, CanonicalMetric}, ::Any...) where {n,k}
-
-function log!(::MetricManifold{ℝ, Stiefel{n,k,ℝ}, CanonicalMetric}, p, q; tolerance=1e-9, maxiter=1e5) where{n,k}
-    M = p'*q
-    QR = qr(q-p*M)
-    V = qr([M;Matrix(QR.R)]).Q*Matrix(I, 2*k,2*k)
-    S = svd(V[(k+1):2*k, (k+1):2k]) #bottom right corner
-    V[:,(k+1):2*k] = V[:,(k+1):2*k]*(S.V*S.U')
+log(::MetricManifold{ℝ,Stiefel{n,k,ℝ},CanonicalMetric}, ::Any...) where {n,k}
+
+function log!(
+    ::MetricManifold{ℝ,Stiefel{n,k,ℝ},CanonicalMetric},
+    p,
+    q;
+    tolerance=1e-9,
+    maxiter=1e5,
+) where {n,k}
+    M = p' * q
+    QR = qr(q - p * M)
+    V = qr([M; Matrix(QR.R)]).Q * Matrix(I, 2 * k, 2 * k)
+    S = svd(V[(k + 1):(2 * k), (k + 1):(2k)]) #bottom right corner
+    V[:, (k + 1):(2 * k)] = V[:, (k + 1):(2 * k)] * (S.V * S.U')
     LV = log(V)
-    C = view(LV, (k+1):2*k, (k+1):2*k)
-    i=0
+    C = view(LV, (k + 1):(2 * k), (k + 1):(2 * k))
+    i = 0
     while (i < maxiter) && (norm(C) > tolerance)
         LV = log(V)
-        V[:, (k+1):2*k] *= exp(-C)
+        V[:, (k + 1):(2 * k)] *= exp(-C)
     end
-    X .= p*LV[1:k,1:k] + QR.Q*LV[(k+1):2*k, 1:k]
+    X .= p * LV[1:k, 1:k] + QR.Q * LV[(k + 1):(2 * k), 1:k]
     return X
 end
diff --git a/src/manifolds/StiefelEuclideanMetric.jl b/src/manifolds/StiefelEuclideanMetric.jl
index c49ba01b88..7c578f7a84 100644
--- a/src/manifolds/StiefelEuclideanMetric.jl
+++ b/src/manifolds/StiefelEuclideanMetric.jl
@@ -1,10 +1,6 @@
 default_metric_dispatch(::Stiefel, ::EuclideanMetric) = Val(true)
 
-function decorator_transparent_dispatch(
-    ::typeof(distance),
-    ::Stiefel,
-    args...,
-)
+function decorator_transparent_dispatch(::typeof(distance), ::Stiefel, args...)
     return Val(:intransparent)
 end
 

From e79a1b5e0cad28b29c84e5160bdea9b4071ac3b8 Mon Sep 17 00:00:00 2001
From: Ronny <git@ronnybergmann.net>
Date: Sat, 6 Mar 2021 18:35:17 +0100
Subject: [PATCH 03/32] undo a commented our part in docs.

---
 docs/make.jl | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/docs/make.jl b/docs/make.jl
index 0319e66897..c4aa2186dd 100644
--- a/docs/make.jl
+++ b/docs/make.jl
@@ -46,7 +46,7 @@ makedocs(
                 "Rotations" => "manifolds/rotations.md",
                 "Skew-symmetric matrices" => "manifolds/skewsymmetric.md",
                 "Spectrahedron" => "manifolds/spectrahedron.md",
-#                "Sphere" => "manifolds/sphere.md",
+                "Sphere" => "manifolds/sphere.md",
                 "Stiefel" => "manifolds/stiefel.md",
                 "Symmetric matrices" => "manifolds/symmetric.md",
                 "Symmetric positive definite" =>

From c993f3ba79ea11fc3b46dae3ca282bf621166adf Mon Sep 17 00:00:00 2001
From: Ronny Bergmann <git@ronnybergmann.net>
Date: Sun, 7 Mar 2021 08:40:02 +0100
Subject: [PATCH 04/32] Apply suggestions from code review

Co-authored-by: Seth Axen <seth.axen@gmail.com>
---
 src/manifolds/StiefelCanonicalMetric.jl | 20 ++++++++++++--------
 src/manifolds/StiefelEuclideanMetric.jl | 14 ++++++--------
 2 files changed, 18 insertions(+), 16 deletions(-)

diff --git a/src/manifolds/StiefelCanonicalMetric.jl b/src/manifolds/StiefelCanonicalMetric.jl
index aa6d1d1d10..63d93b8ba5 100644
--- a/src/manifolds/StiefelCanonicalMetric.jl
+++ b/src/manifolds/StiefelCanonicalMetric.jl
@@ -54,8 +54,8 @@ exp(::MetricManifold{ℝ,Stiefel{n,k,ℝ},CanonicalMetric}, ::Any...) where {n,k
 function exp!(::MetricManifold{ℝ,Stiefel{n,k,ℝ},CanonicalMetric}, Q, p, X) where {n,k}
     A = p' * X
     QR = qr(X - p * A)
-    BC_ext = exp([A -QR.R'; QR.R, zeros(k, k)])
-    q .= [p Matrix(QR.Q)] * BC_ext[:, 1:k]
+    BC_ext = exp([A -QR.R'; QR.R 0*I])
+    q .= [p Matrix(QR.Q)] * @view(BC_ext[:, 1:k])
     return q
 end
 
@@ -70,7 +70,7 @@ g_p(X,Y) = \operatorname{tr}\bigl( X^{\mathrm{T}}(I_n - \frac{1}{2}pp^{\mathrm{T
 ```
 """
 function inner(::MetricManifold{ℝ,Stiefel{n,k,ℝ},X,CanonicalMetric}, p, X, Y) where {n,k}
-    return tr(X' * (Matrix(I, n, n) - 0.5 * p * p') * Y)
+    return dot(X, (Y .- (p * (p'Y)) ./ 2))
 end
 
 @doc raw"""
@@ -93,6 +93,7 @@ log(::MetricManifold{ℝ,Stiefel{n,k,ℝ},CanonicalMetric}, ::Any...) where {n,k
 
 function log!(
     ::MetricManifold{ℝ,Stiefel{n,k,ℝ},CanonicalMetric},
+    X,
     p,
     q;
     tolerance=1e-9,
@@ -100,16 +101,19 @@ function log!(
 ) where {n,k}
     M = p' * q
     QR = qr(q - p * M)
-    V = qr([M; Matrix(QR.R)]).Q * Matrix(I, 2 * k, 2 * k)
-    S = svd(V[(k + 1):(2 * k), (k + 1):(2k)]) #bottom right corner
-    V[:, (k + 1):(2 * k)] = V[:, (k + 1):(2 * k)] * (S.V * S.U')
+    V = Matrix(qr([M; QR.R]).Q)
+    Vcorner = @view V[(k + 1):(2 * k), (k + 1):(2k)] #bottom right corner
+    S = svd!(Vcorner)
+    mul!(Vcorner, Vcorner * S.V, S.U')
     LV = log(V)
     C = view(LV, (k + 1):(2 * k), (k + 1):(2 * k))
+    expnC = exp(-C)
     i = 0
     while (i < maxiter) && (norm(C) > tolerance)
         LV = log(V)
-        V[:, (k + 1):(2 * k)] *= exp(-C)
+        copyto!(Vcorner, Vcorner*expnC)
     end
-    X .= p * LV[1:k, 1:k] + QR.Q * LV[(k + 1):(2 * k), 1:k]
+    mul!(X, p, @view(LV[1:k, 1:k]))
+    mul!(X, QR.Q, @view(LV[(k + 1):(2 * k), 1:k]), true, true)
     return X
 end
diff --git a/src/manifolds/StiefelEuclideanMetric.jl b/src/manifolds/StiefelEuclideanMetric.jl
index 7c578f7a84..95968e888d 100644
--- a/src/manifolds/StiefelEuclideanMetric.jl
+++ b/src/manifolds/StiefelEuclideanMetric.jl
@@ -29,12 +29,8 @@ $0_k$ are the identity matrix and the zero matrix of dimension $k × k$, respect
 exp(::Stiefel, ::Any...)
 
 function exp!(::Stiefel{n,k}, q, p, X) where {n,k}
-    return copyto!(
-        q,
-        [p X] *
-        exp([p'X -X'*X; one(zeros(eltype(p), k, k)) p'*X]) *
-        [exp(-p'X); zeros(eltype(p), k, k)],
-    )
+    A = p'*X
+    return copyto!(q, [p X] * exp([A -X'*X; I p'*X]) * [exp(-A); 0*I])
 end
 
 @doc raw"""
@@ -74,7 +70,7 @@ function get_coordinates!(
     B::DefaultOrthonormalBasis{ℝ},
 ) where {n,k}
     V = get_vectors(M, p, B)
-    c .= [inner(M, p, v, X) for v in V]
+    c .= inner.(Ref(M), Ref(p), V, Ref(X))
     return c
 end
 
@@ -135,7 +131,9 @@ project(::Stiefel, ::Any...)
 
 function project!(::Stiefel, Y, p, X)
     A = p' * X
-    copyto!(Y, X - p * Hermitian((A + A') / 2))
+    T = eltype(Y)
+    copyto!(Y, X)
+    mul!(Y, p, A + A', T(-0.5), true)
     return Y
 end
 

From f36ab3ef29e3a2af8f972ffb0b495e0637b61719 Mon Sep 17 00:00:00 2001
From: Ronny Bergmann <git@ronnybergmann.net>
Date: Mon, 8 Mar 2021 07:28:38 +0100
Subject: [PATCH 05/32] Apply suggestions from code review

Co-authored-by: Seth Axen <seth.axen@gmail.com>
---
 src/manifolds/StiefelCanonicalMetric.jl | 8 +++++++-
 src/manifolds/StiefelEuclideanMetric.jl | 2 +-
 2 files changed, 8 insertions(+), 2 deletions(-)

diff --git a/src/manifolds/StiefelCanonicalMetric.jl b/src/manifolds/StiefelCanonicalMetric.jl
index 63d93b8ba5..16554ce0b7 100644
--- a/src/manifolds/StiefelCanonicalMetric.jl
+++ b/src/manifolds/StiefelCanonicalMetric.jl
@@ -53,6 +53,7 @@ exp(::MetricManifold{ℝ,Stiefel{n,k,ℝ},CanonicalMetric}, ::Any...) where {n,k
 
 function exp!(::MetricManifold{ℝ,Stiefel{n,k,ℝ},CanonicalMetric}, Q, p, X) where {n,k}
     A = p' * X
+    n == k && return mul!(q, p, exp(A))
     QR = qr(X - p * A)
     BC_ext = exp([A -QR.R'; QR.R 0*I])
     q .= [p Matrix(QR.Q)] * @view(BC_ext[:, 1:k])
@@ -70,7 +71,12 @@ g_p(X,Y) = \operatorname{tr}\bigl( X^{\mathrm{T}}(I_n - \frac{1}{2}pp^{\mathrm{T
 ```
 """
 function inner(::MetricManifold{ℝ,Stiefel{n,k,ℝ},X,CanonicalMetric}, p, X, Y) where {n,k}
-    return dot(X, (Y .- (p * (p'Y)) ./ 2))
+    T = Base.promote_eltype(p, X, Y)
+    s = T(dot(X, Y))
+    if n != k
+        s -= dot(p'X, p'Y) / 2
+    end
+    return s
 end
 
 @doc raw"""
diff --git a/src/manifolds/StiefelEuclideanMetric.jl b/src/manifolds/StiefelEuclideanMetric.jl
index 95968e888d..f093c37e9a 100644
--- a/src/manifolds/StiefelEuclideanMetric.jl
+++ b/src/manifolds/StiefelEuclideanMetric.jl
@@ -30,7 +30,7 @@ exp(::Stiefel, ::Any...)
 
 function exp!(::Stiefel{n,k}, q, p, X) where {n,k}
     A = p'*X
-    return copyto!(q, [p X] * exp([A -X'*X; I p'*X]) * [exp(-A); 0*I])
+    return copyto!(q, [p X] * exp([A -X'*X; I A]) * [exp(-A); 0*I])
 end
 
 @doc raw"""

From edc96c95f01ba196eada20ad2476dcbc6ec0ae8c Mon Sep 17 00:00:00 2001
From: Ronny <git@ronnybergmann.net>
Date: Mon, 8 Mar 2021 07:45:54 +0100
Subject: [PATCH 06/32] runs formatter.

---
 src/manifolds/StiefelCanonicalMetric.jl | 2 +-
 src/manifolds/StiefelEuclideanMetric.jl | 4 ++--
 2 files changed, 3 insertions(+), 3 deletions(-)

diff --git a/src/manifolds/StiefelCanonicalMetric.jl b/src/manifolds/StiefelCanonicalMetric.jl
index 16554ce0b7..14aa20fded 100644
--- a/src/manifolds/StiefelCanonicalMetric.jl
+++ b/src/manifolds/StiefelCanonicalMetric.jl
@@ -117,7 +117,7 @@ function log!(
     i = 0
     while (i < maxiter) && (norm(C) > tolerance)
         LV = log(V)
-        copyto!(Vcorner, Vcorner*expnC)
+        copyto!(Vcorner, Vcorner * expnC)
     end
     mul!(X, p, @view(LV[1:k, 1:k]))
     mul!(X, QR.Q, @view(LV[(k + 1):(2 * k), 1:k]), true, true)
diff --git a/src/manifolds/StiefelEuclideanMetric.jl b/src/manifolds/StiefelEuclideanMetric.jl
index f093c37e9a..d1d359e7e8 100644
--- a/src/manifolds/StiefelEuclideanMetric.jl
+++ b/src/manifolds/StiefelEuclideanMetric.jl
@@ -29,8 +29,8 @@ $0_k$ are the identity matrix and the zero matrix of dimension $k × k$, respect
 exp(::Stiefel, ::Any...)
 
 function exp!(::Stiefel{n,k}, q, p, X) where {n,k}
-    A = p'*X
-    return copyto!(q, [p X] * exp([A -X'*X; I A]) * [exp(-A); 0*I])
+    A = p' * X
+    return copyto!(q, [p X] * exp([A -X'*X; I A]) * [exp(-A); 0 * I])
 end
 
 @doc raw"""

From e7b914ef3af99b6eb03cebf7261591647a468fd2 Mon Sep 17 00:00:00 2001
From: Ronny <git@ronnybergmann.net>
Date: Mon, 8 Mar 2021 07:54:58 +0100
Subject: [PATCH 07/32] Fix plots (for documentation and testing) to 1.10.5,
 since it introduces a bug (see
 https://github.com/JuliaPlots/Plots.jl/issues/3341) in the next version

---
 Project.toml | 1 +
 1 file changed, 1 insertion(+)

diff --git a/Project.toml b/Project.toml
index dc8e68d76b..19d5d8bd91 100644
--- a/Project.toml
+++ b/Project.toml
@@ -31,6 +31,7 @@ FiniteDifferences = "0.12"
 HybridArrays = "0.4"
 LightGraphs = "1"
 ManifoldsBase = "0.10.1"
+Plots = "= 1.10.5"
 RecipesBase = "1.1"
 Requires = "0.5, 1"
 SimpleWeightedGraphs = "1"

From 3950004f35c647e4efc6198b695a00579cc889bb Mon Sep 17 00:00:00 2001
From: Ronny <git@ronnybergmann.net>
Date: Mon, 8 Mar 2021 08:08:40 +0100
Subject: [PATCH 08/32] Losen constraints for Julia 1.4

---
 Project.toml | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/Project.toml b/Project.toml
index 19d5d8bd91..8f915b6031 100644
--- a/Project.toml
+++ b/Project.toml
@@ -31,7 +31,7 @@ FiniteDifferences = "0.12"
 HybridArrays = "0.4"
 LightGraphs = "1"
 ManifoldsBase = "0.10.1"
-Plots = "= 1.10.5"
+Plots = "1.6, = 1.10.5"
 RecipesBase = "1.1"
 Requires = "0.5, 1"
 SimpleWeightedGraphs = "1"

From d9008f31491dd2a4e9cdc38714f873154ffb8d02 Mon Sep 17 00:00:00 2001
From: Ronny <git@ronnybergmann.net>
Date: Mon, 8 Mar 2021 08:13:17 +0100
Subject: [PATCH 09/32] fix latex equation.

---
 docs/src/manifolds/stiefel.md | 5 +++--
 1 file changed, 3 insertions(+), 2 deletions(-)

diff --git a/docs/src/manifolds/stiefel.md b/docs/src/manifolds/stiefel.md
index e07d10d958..a23ef419a5 100644
--- a/docs/src/manifolds/stiefel.md
+++ b/docs/src/manifolds/stiefel.md
@@ -31,8 +31,9 @@ Any ``X∈T_p\mathcal M``, ``p∈\mathcal M``, can be written as
 ```math
 X = pA + (I_n-pp^{\mathrm{T}})B,
 \quad
-A ∈ ℝ^{p\imes p} \text{skew-symmetric},
-B ∈ ℝ^{n\times p} \text{arbitrary.}
+A ∈ ℝ^{p×p} \text{ skew-symmetric},
+\quad
+B ∈ ℝ^{n×p} \text{ arbitrary.}
 ```
 
 In the Euclidean case, the elements from ``A`` are counted twice (i.e. weighted with a factor of 2).

From 9c61049ebdcb756a55d530af6beafc9a4d7ce5cb Mon Sep 17 00:00:00 2001
From: Ronny <git@ronnybergmann.net>
Date: Mon, 8 Mar 2021 10:19:45 +0100
Subject: [PATCH 10/32] Fix documentation.

---
 docs/Project.toml                       |  2 +-
 src/Manifolds.jl                        |  1 +
 src/manifolds/StiefelCanonicalMetric.jl | 10 +++++-----
 3 files changed, 7 insertions(+), 6 deletions(-)

diff --git a/docs/Project.toml b/docs/Project.toml
index 617ae69605..54cfbeed2f 100644
--- a/docs/Project.toml
+++ b/docs/Project.toml
@@ -17,6 +17,6 @@ Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40"
 Documenter = "0.24, 0.25"
 HybridArrays = "0.4"
 ManifoldsBase = "0.10"
-Plots = "1.7"
+Plots = "1.6, 1.7, = 1.10.5"
 StaticArrays = "1.0"
 PyPlot = "2.9"
diff --git a/src/Manifolds.jl b/src/Manifolds.jl
index e29788f356..0fa4bc6d36 100644
--- a/src/Manifolds.jl
+++ b/src/Manifolds.jl
@@ -161,6 +161,7 @@ include("manifolds/SkewSymmetric.jl")
 include("manifolds/Spectrahedron.jl")
 include("manifolds/Stiefel.jl")
 include("manifolds/StiefelEuclideanMetric.jl")
+include("manifolds/StiefelCanonicalMetric.jl")
 include("manifolds/Sphere.jl")
 include("manifolds/SphereSymmetricMatrices.jl")
 include("manifolds/Symmetric.jl")
diff --git a/src/manifolds/StiefelCanonicalMetric.jl b/src/manifolds/StiefelCanonicalMetric.jl
index 14aa20fded..b9e42c2e2d 100644
--- a/src/manifolds/StiefelCanonicalMetric.jl
+++ b/src/manifolds/StiefelCanonicalMetric.jl
@@ -26,8 +26,8 @@ respect to ``p``, i.e.
 X = pp^{\mathrm{T}}X + (I_n-pp^{\mathrm{T}})X,
 ```
 where ``I_n`` is the ``n\times n`` identity matrix.
-We introduce ``A=p^{\mathrm{T}}X`` and ``QR = (I_n-pp^{\mathrm{T}})X``` the qr decomposition
-of the vertical component. Then using the matrix exponential ``operatorname{Exp}`` we introduce ``B`` and ``C`` as
+We introduce ``A=p^{\mathrm{T}}X`` and ``QR = (I_n-pp^{\mathrm{T}})X`` the `qr` decomposition
+of the vertical component. Then using the matrix exponential ``\operatorname{Exp}`` we introduce ``B`` and ``C`` as
 
 ```math
 \begin{pmatrix}
@@ -66,11 +66,11 @@ end
 compute the inner procuct on the [`Stiefel`](@ref) manifold with respect to the
 [`CanonicalMetric`](@ref). The formula reads
 
-```
+```math
 g_p(X,Y) = \operatorname{tr}\bigl( X^{\mathrm{T}}(I_n - \frac{1}{2}pp^{\mathrm{T}})Y \bigr).
 ```
 """
-function inner(::MetricManifold{ℝ,Stiefel{n,k,ℝ},X,CanonicalMetric}, p, X, Y) where {n,k}
+function inner(::MetricManifold{ℝ,Stiefel{n,k,ℝ},CanonicalMetric}, p, X, Y) where {n,k}
     T = Base.promote_eltype(p, X, Y)
     s = T(dot(X, Y))
     if n != k
@@ -85,7 +85,7 @@ end
 
 Compute the logarithmic map on the [`Stiefel`](@ref)`(n,k)` manifold with respect to the [`CanonicalMetric`](@ref)
 using a matrix-algebraic based approach to an iterative inversion of the formula of the
-[`exp`](@ref exp(::MetricManifold{ℝ, Stiefel{n,k,ℝ}, CanonicalMetric}, ::Any...)).
+[`exp`](@ref exp(::MetricManifold{ℝ, Stiefel{n,k,ℝ}, CanonicalMetric}, ::Any...) where {n,k}).
 
 The algorithm is derived in[^Zimmermann2017].
 

From 0a8aae0bd8a73c85e24fa85c81180ed7080ed217 Mon Sep 17 00:00:00 2001
From: Ronny <git@ronnybergmann.net>
Date: Mon, 8 Mar 2021 10:36:34 +0100
Subject: [PATCH 11/32] changing compats again, hopefully this works?

---
 Project.toml      | 2 +-
 docs/Project.toml | 2 +-
 2 files changed, 2 insertions(+), 2 deletions(-)

diff --git a/Project.toml b/Project.toml
index 8f915b6031..de11bf48ef 100644
--- a/Project.toml
+++ b/Project.toml
@@ -31,7 +31,7 @@ FiniteDifferences = "0.12"
 HybridArrays = "0.4"
 LightGraphs = "1"
 ManifoldsBase = "0.10.1"
-Plots = "1.6, = 1.10.5"
+Plots = "~1.6, =1.10.5"
 RecipesBase = "1.1"
 Requires = "0.5, 1"
 SimpleWeightedGraphs = "1"
diff --git a/docs/Project.toml b/docs/Project.toml
index 54cfbeed2f..ba550568d5 100644
--- a/docs/Project.toml
+++ b/docs/Project.toml
@@ -17,6 +17,6 @@ Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40"
 Documenter = "0.24, 0.25"
 HybridArrays = "0.4"
 ManifoldsBase = "0.10"
-Plots = "1.6, 1.7, = 1.10.5"
+Plots = "= 1.10.5"
 StaticArrays = "1.0"
 PyPlot = "2.9"

From f97f37fa852379ef81860a0b2ca2d6d249676fbd Mon Sep 17 00:00:00 2001
From: Ronny <git@ronnybergmann.net>
Date: Mon, 8 Mar 2021 14:04:10 +0100
Subject: [PATCH 12/32] starts testing.

---
 src/Manifolds.jl                        |  1 +
 src/manifolds/StiefelCanonicalMetric.jl |  9 +++++----
 test/stiefel.jl                         | 18 ++++++++++++++++++
 3 files changed, 24 insertions(+), 4 deletions(-)

diff --git a/src/Manifolds.jl b/src/Manifolds.jl
index 0fa4bc6d36..2bfbd1ad9b 100644
--- a/src/Manifolds.jl
+++ b/src/Manifolds.jl
@@ -333,6 +333,7 @@ export Metric,
     MinkowskiMetric,
     PowerMetric,
     ProductMetric,
+    CanonicalMetric,
     MetricManifold
 export AbstractEmbeddingType, AbstractIsometricEmbeddingType
 export DefaultEmbeddingType, DefaultIsometricEmbeddingType, TransparentIsometricEmbedding
diff --git a/src/manifolds/StiefelCanonicalMetric.jl b/src/manifolds/StiefelCanonicalMetric.jl
index b9e42c2e2d..0e0d8755be 100644
--- a/src/manifolds/StiefelCanonicalMetric.jl
+++ b/src/manifolds/StiefelCanonicalMetric.jl
@@ -51,7 +51,7 @@ For more details, see [^EdelmanAriasSmith1998][^Zimmermann2017].
 """
 exp(::MetricManifold{ℝ,Stiefel{n,k,ℝ},CanonicalMetric}, ::Any...) where {n,k}
 
-function exp!(::MetricManifold{ℝ,Stiefel{n,k,ℝ},CanonicalMetric}, Q, p, X) where {n,k}
+function exp!(::MetricManifold{ℝ,Stiefel{n,k,ℝ},CanonicalMetric}, q, p, X) where {n,k}
     A = p' * X
     n == k && return mul!(q, p, exp(A))
     QR = qr(X - p * A)
@@ -107,11 +107,12 @@ function log!(
 ) where {n,k}
     M = p' * q
     QR = qr(q - p * M)
-    V = Matrix(qr([M; QR.R]).Q)
+    V = Matrix(qr([M; QR.R]).Q*I)
     Vcorner = @view V[(k + 1):(2 * k), (k + 1):(2k)] #bottom right corner
     S = svd!(Vcorner)
-    mul!(Vcorner, Vcorner * S.V, S.U')
-    LV = log(V)
+    mul!(Vcorner, Vcorner * S.U, S.V')
+    V[:,1:k] .= [M; QR.R]
+    LV = real(log(V))
     C = view(LV, (k + 1):(2 * k), (k + 1):(2 * k))
     expnC = exp(-C)
     i = 0
diff --git a/test/stiefel.jl b/test/stiefel.jl
index 561830cc44..7669a2180b 100644
--- a/test/stiefel.jl
+++ b/test/stiefel.jl
@@ -1,11 +1,15 @@
 include("utils.jl")
 
+using Manifolds: default_metric_dispatch
+
 @testset "Stiefel" begin
     @testset "Real" begin
         M = Stiefel(3, 2)
+        M2 = MetricManifold(M, EuclideanMetric())
         @testset "Basics" begin
             @test repr(M) == "Stiefel(3, 2, ℝ)"
             x = [1.0 0.0; 0.0 1.0; 0.0 0.0]
+            @test (@inferred default_metric_dispatch(M2)) === Val(true)
             @test representation_size(M) == (3, 2)
             @test manifold_dimension(M) == 3
             base_manifold(M) === M
@@ -22,6 +26,7 @@ include("utils.jl")
                 1 * im * zero_tangent_vector(M, x),
                 true,
             )
+
         end
         @testset "Embedding and Projection" begin
             x = [1.0 0.0; 0.0 1.0; 0.0 0.0]
@@ -92,6 +97,7 @@ include("utils.jl")
             x = [1.0 0.0; 0.0 1.0; 0.0 0.0]
             y = exp(M, x, [0.0 0.0; 0.0 0.0; 1.0 1.0])
             z = exp(M, x, [0.0 0.0; 0.0 0.0; -1.0 1.0])
+            @test_throws ErrorException distance(M, x, y)
             @test isapprox(
                 M,
                 retract(
@@ -296,4 +302,16 @@ include("utils.jl")
         q2 = retract(M, p, X, r2)
         @test is_manifold_point(M, q2)
     end
+
+    @testset "Canonical Metric" begin
+        M3 = MetricManifold(Stiefel(3,2), CanonicalMetric())
+        p = [1.0 0.0; 0.0 1.0; 0.0 0.0]
+        X = [0.0 0.0; 0.0 0.0; 1.0 1.0]
+        q = exp(M3, p, X)
+        Y = [0.0 0.0; 0.0 0.0; -1.0 1.0]
+        r = exp(M3, p, Y)
+        @test isapprox(M3, p, log(M3, p, q), X)
+        @test isapprox(M3, p, log(M3, p, r), Y)
+        @test inner(M3, p, X, Y) == 0
+    end
 end

From cee6df24075641798a0a73c99b723145c2502120 Mon Sep 17 00:00:00 2001
From: Ronny <git@ronnybergmann.net>
Date: Mon, 8 Mar 2021 14:11:25 +0100
Subject: [PATCH 13/32] adopts one view, such that the multiplication is well
 defined. Finishes testing.

---
 src/manifolds/StiefelCanonicalMetric.jl | 7 ++++---
 1 file changed, 4 insertions(+), 3 deletions(-)

diff --git a/src/manifolds/StiefelCanonicalMetric.jl b/src/manifolds/StiefelCanonicalMetric.jl
index 0e0d8755be..c5044a7fb6 100644
--- a/src/manifolds/StiefelCanonicalMetric.jl
+++ b/src/manifolds/StiefelCanonicalMetric.jl
@@ -107,11 +107,11 @@ function log!(
 ) where {n,k}
     M = p' * q
     QR = qr(q - p * M)
-    V = Matrix(qr([M; QR.R]).Q*I)
+    V = Matrix(qr([M; QR.R]).Q * I)
     Vcorner = @view V[(k + 1):(2 * k), (k + 1):(2k)] #bottom right corner
     S = svd!(Vcorner)
     mul!(Vcorner, Vcorner * S.U, S.V')
-    V[:,1:k] .= [M; QR.R]
+    V[:, 1:k] .= [M; QR.R]
     LV = real(log(V))
     C = view(LV, (k + 1):(2 * k), (k + 1):(2 * k))
     expnC = exp(-C)
@@ -121,6 +121,7 @@ function log!(
         copyto!(Vcorner, Vcorner * expnC)
     end
     mul!(X, p, @view(LV[1:k, 1:k]))
-    mul!(X, QR.Q, @view(LV[(k + 1):(2 * k), 1:k]), true, true)
+    # force the first - Q - to be the reduced form, not the full matrix
+    mul!(X, @view(QR.Q[:, 1:k]), @view(LV[(k + 1):(2 * k), 1:k]), true, true)
     return X
 end

From 018c288fd0bc7be001b853f6b7efa3b86b652a54 Mon Sep 17 00:00:00 2001
From: Ronny <git@ronnybergmann.net>
Date: Mon, 8 Mar 2021 15:02:20 +0100
Subject: [PATCH 14/32] Fix spacing.

---
 test/stiefel.jl | 3 +--
 1 file changed, 1 insertion(+), 2 deletions(-)

diff --git a/test/stiefel.jl b/test/stiefel.jl
index 7669a2180b..030c040483 100644
--- a/test/stiefel.jl
+++ b/test/stiefel.jl
@@ -26,7 +26,6 @@ using Manifolds: default_metric_dispatch
                 1 * im * zero_tangent_vector(M, x),
                 true,
             )
-
         end
         @testset "Embedding and Projection" begin
             x = [1.0 0.0; 0.0 1.0; 0.0 0.0]
@@ -304,7 +303,7 @@ using Manifolds: default_metric_dispatch
     end
 
     @testset "Canonical Metric" begin
-        M3 = MetricManifold(Stiefel(3,2), CanonicalMetric())
+        M3 = MetricManifold(Stiefel(3, 2), CanonicalMetric())
         p = [1.0 0.0; 0.0 1.0; 0.0 0.0]
         X = [0.0 0.0; 0.0 0.0; 1.0 1.0]
         q = exp(M3, p, X)

From c89c7636ea47a256276944f2c05f23bab839e36b Mon Sep 17 00:00:00 2001
From: Ronny <git@ronnybergmann.net>
Date: Mon, 8 Mar 2021 17:12:44 +0100
Subject: [PATCH 15/32] fixes a few bugs, starts extended testing that includes
 the iterations.

---
 src/manifolds/StiefelCanonicalMetric.jl | 12 ++++++++----
 test/stiefel.jl                         | 10 +++++++++-
 2 files changed, 17 insertions(+), 5 deletions(-)

diff --git a/src/manifolds/StiefelCanonicalMetric.jl b/src/manifolds/StiefelCanonicalMetric.jl
index c5044a7fb6..d0432d0741 100644
--- a/src/manifolds/StiefelCanonicalMetric.jl
+++ b/src/manifolds/StiefelCanonicalMetric.jl
@@ -107,19 +107,23 @@ function log!(
 ) where {n,k}
     M = p' * q
     QR = qr(q - p * M)
-    V = Matrix(qr([M; QR.R]).Q * I)
+    V = Complex.(Matrix(qr([M; QR.R]).Q * I))
     Vcorner = @view V[(k + 1):(2 * k), (k + 1):(2k)] #bottom right corner
+    Vpcols = @view V[:, (k + 1):(2 * k)]
     S = svd!(Vcorner)
     mul!(Vcorner, Vcorner * S.U, S.V')
     V[:, 1:k] .= [M; QR.R]
-    LV = real(log(V))
+    LV = log(V)
     C = view(LV, (k + 1):(2 * k), (k + 1):(2 * k))
     expnC = exp(-C)
     i = 0
     while (i < maxiter) && (norm(C) > tolerance)
-        LV = log(V)
-        copyto!(Vcorner, Vcorner * expnC)
+        i = i + 1
+        LV .= log(V)
+        expnC .= exp(-C)
+        copyto!(Vpcols, Vpcols * expnC)
     end
+    LV .= real.(LV)
     mul!(X, p, @view(LV[1:k, 1:k]))
     # force the first - Q - to be the reduced form, not the full matrix
     mul!(X, @view(QR.Q[:, 1:k]), @view(LV[(k + 1):(2 * k), 1:k]), true, true)
diff --git a/test/stiefel.jl b/test/stiefel.jl
index 030c040483..85a29a4f3e 100644
--- a/test/stiefel.jl
+++ b/test/stiefel.jl
@@ -311,6 +311,14 @@ using Manifolds: default_metric_dispatch
         r = exp(M3, p, Y)
         @test isapprox(M3, p, log(M3, p, q), X)
         @test isapprox(M3, p, log(M3, p, r), Y)
+        Z = [0.0 sqrt(2); -sqrt(2) 0.0; -sqrt(2) sqrt(2)]
+        s = exp(M3, p, Z)
+        @test is_tangent_vector(M3, p, s)
         @test inner(M3, p, X, Y) == 0
+        a = 1.4
+        Z = [0.0 2a; -2a 0.0; -a/2 a/2]
+        s = exp(M3, p, Z)
+        Z2 = log(M3, p, s)
+        @test_broken isapprox(M3, p, Z, Z2)
     end
-end
+end
\ No newline at end of file

From b062c87c61505a8620d103b08fea76ea1900d177 Mon Sep 17 00:00:00 2001
From: Ronny <git@ronnybergmann.net>
Date: Mon, 8 Mar 2021 17:28:45 +0100
Subject: [PATCH 16/32] runs formatter.

---
 test/stiefel.jl | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/test/stiefel.jl b/test/stiefel.jl
index 85a29a4f3e..1bca16239a 100644
--- a/test/stiefel.jl
+++ b/test/stiefel.jl
@@ -321,4 +321,4 @@ using Manifolds: default_metric_dispatch
         Z2 = log(M3, p, s)
         @test_broken isapprox(M3, p, Z, Z2)
     end
-end
\ No newline at end of file
+end

From 50f83cabb9ee036f454884534ca0b72547ca21dd Mon Sep 17 00:00:00 2001
From: Ronny <git@ronnybergmann.net>
Date: Mon, 8 Mar 2021 18:11:21 +0100
Subject: [PATCH 17/32] removes a little bit of unnecessary lines.

---
 test/stiefel.jl | 3 ---
 1 file changed, 3 deletions(-)

diff --git a/test/stiefel.jl b/test/stiefel.jl
index 1bca16239a..4897209768 100644
--- a/test/stiefel.jl
+++ b/test/stiefel.jl
@@ -311,9 +311,6 @@ using Manifolds: default_metric_dispatch
         r = exp(M3, p, Y)
         @test isapprox(M3, p, log(M3, p, q), X)
         @test isapprox(M3, p, log(M3, p, r), Y)
-        Z = [0.0 sqrt(2); -sqrt(2) 0.0; -sqrt(2) sqrt(2)]
-        s = exp(M3, p, Z)
-        @test is_tangent_vector(M3, p, s)
         @test inner(M3, p, X, Y) == 0
         a = 1.4
         Z = [0.0 2a; -2a 0.0; -a/2 a/2]

From 4d359b5f4fcbd9855f32f2a342d8d01fd52ed67b Mon Sep 17 00:00:00 2001
From: Ronny <git@ronnybergmann.net>
Date: Thu, 11 Mar 2021 11:58:18 +0100
Subject: [PATCH 18/32] enforce real part of log and avoid svd! to not change
 V.

---
 src/manifolds/StiefelCanonicalMetric.jl | 13 +++++++------
 test/stiefel.jl                         |  2 +-
 2 files changed, 8 insertions(+), 7 deletions(-)

diff --git a/src/manifolds/StiefelCanonicalMetric.jl b/src/manifolds/StiefelCanonicalMetric.jl
index d0432d0741..71eb0f5095 100644
--- a/src/manifolds/StiefelCanonicalMetric.jl
+++ b/src/manifolds/StiefelCanonicalMetric.jl
@@ -107,19 +107,20 @@ function log!(
 ) where {n,k}
     M = p' * q
     QR = qr(q - p * M)
-    V = Complex.(Matrix(qr([M; QR.R]).Q * I))
+    V = Matrix(qr([M; QR.R]).Q * I)
     Vcorner = @view V[(k + 1):(2 * k), (k + 1):(2k)] #bottom right corner
-    Vpcols = @view V[:, (k + 1):(2 * k)]
-    S = svd!(Vcorner)
-    mul!(Vcorner, Vcorner * S.U, S.V')
+    Vpcols = @view V[:, (k + 1):(2 * k)] #second half of the columns
+    S = svd(Vcorner) # preprocessing: Procrustes
+    mul!(Vpcols, Vpcols * S.U, S.V')
     V[:, 1:k] .= [M; QR.R]
-    LV = log(V)
+
+    LV = real.(log(V)) # this can be replaced by log_safe.
     C = view(LV, (k + 1):(2 * k), (k + 1):(2 * k))
     expnC = exp(-C)
     i = 0
     while (i < maxiter) && (norm(C) > tolerance)
         i = i + 1
-        LV .= log(V)
+        LV .= real.(log(V))
         expnC .= exp(-C)
         copyto!(Vpcols, Vpcols * expnC)
     end
diff --git a/test/stiefel.jl b/test/stiefel.jl
index 4897209768..ed82867de9 100644
--- a/test/stiefel.jl
+++ b/test/stiefel.jl
@@ -316,6 +316,6 @@ using Manifolds: default_metric_dispatch
         Z = [0.0 2a; -2a 0.0; -a/2 a/2]
         s = exp(M3, p, Z)
         Z2 = log(M3, p, s)
-        @test_broken isapprox(M3, p, Z, Z2)
+        @test isapprox(M3, p, Z, Z2)
     end
 end

From 7103ada6220dcd2299935ff349a34e2940ba74f9 Mon Sep 17 00:00:00 2001
From: Ronny Bergmann <git@ronnybergmann.net>
Date: Thu, 11 Mar 2021 19:23:46 +0100
Subject: [PATCH 19/32] Update src/manifolds/StiefelCanonicalMetric.jl

Co-authored-by: Mateusz Baran <mateuszbaran89@gmail.com>
---
 src/manifolds/StiefelCanonicalMetric.jl | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/src/manifolds/StiefelCanonicalMetric.jl b/src/manifolds/StiefelCanonicalMetric.jl
index 71eb0f5095..e48db1d2db 100644
--- a/src/manifolds/StiefelCanonicalMetric.jl
+++ b/src/manifolds/StiefelCanonicalMetric.jl
@@ -61,7 +61,7 @@ function exp!(::MetricManifold{ℝ,Stiefel{n,k,ℝ},CanonicalMetric}, q, p, X) w
 end
 
 @doc raw"""
-    inner(M::MetricManifold{ℝ, Stiefel{n,k,ℝ}, X, CanonicalMetric}, p,, X, Y)
+    inner(M::MetricManifold{ℝ, Stiefel{n,k,ℝ}, X, CanonicalMetric}, p, X, Y)
 
 compute the inner procuct on the [`Stiefel`](@ref) manifold with respect to the
 [`CanonicalMetric`](@ref). The formula reads

From 102018bf6efb3185aecfd3b1bcc3c6aa14929852 Mon Sep 17 00:00:00 2001
From: Ronny <git@ronnybergmann.net>
Date: Thu, 11 Mar 2021 20:59:11 +0100
Subject: [PATCH 20/32] introduce distance and kwargs for the log, extend
 tests.

---
 src/Manifolds.jl                        |  1 +
 src/manifolds/StiefelCanonicalMetric.jl | 29 ++++++++++++++++++++++---
 test/stiefel.jl                         |  2 ++
 3 files changed, 29 insertions(+), 3 deletions(-)

diff --git a/src/Manifolds.jl b/src/Manifolds.jl
index 2bfbd1ad9b..a19e9a10a7 100644
--- a/src/Manifolds.jl
+++ b/src/Manifolds.jl
@@ -40,6 +40,7 @@ import ManifoldsBase:
     is_tangent_vector,
     inverse_retract,
     inverse_retract!,
+    log, #for extension, e.g. in Stiefel.
     log!,
     manifold_dimension,
     mid_point,
diff --git a/src/manifolds/StiefelCanonicalMetric.jl b/src/manifolds/StiefelCanonicalMetric.jl
index 71eb0f5095..46a2b6221c 100644
--- a/src/manifolds/StiefelCanonicalMetric.jl
+++ b/src/manifolds/StiefelCanonicalMetric.jl
@@ -13,6 +13,10 @@ manifold, see for example[^EdelmanAriasSmith1998].
 """
 struct CanonicalMetric <: RiemannianMetric end
 
+function distance(M::MetricManifold{ℝ,Stiefel{n,k,ℝ},CanonicalMetric}, q, p) where {n,k}
+    return norm(M, p, log(M, p, q))
+end
+
 @doc raw"""
     q = exp(M::MetricManifold{ℝ, Stiefel{n,k,ℝ}, CanonicalMetric}, p, X)
     exp!(M::MetricManifold{ℝ, Stiefel{n,k,ℝ}, q, CanonicalMetric}, p, X)
@@ -80,8 +84,8 @@ function inner(::MetricManifold{ℝ,Stiefel{n,k,ℝ},CanonicalMetric}, p, X, Y)
 end
 
 @doc raw"""
-    X = log(M::MetricManifold{ℝ, Stiefel{n,k,ℝ}, CanonicalMetric}, p, q)
-    log!(M::MetricManifold{ℝ, Stiefel{n,k,ℝ}, X, CanonicalMetric}, p, q)
+    X = log(M::MetricManifold{ℝ, Stiefel{n,k,ℝ}, CanonicalMetric}, p, q; kwargs..)
+    log!(M::MetricManifold{ℝ, Stiefel{n,k,ℝ}, X, CanonicalMetric}, p, q; kwargs...)
 
 Compute the logarithmic map on the [`Stiefel`](@ref)`(n,k)` manifold with respect to the [`CanonicalMetric`](@ref)
 using a matrix-algebraic based approach to an iterative inversion of the formula of the
@@ -89,6 +93,11 @@ using a matrix-algebraic based approach to an iterative inversion of the formula
 
 The algorithm is derived in[^Zimmermann2017].
 
+# Keyword arguments
+
+Both `maxiter`(`=1e5`) and `tolerance`(`=1e-9`) can be given to the iterative process as
+stoppping criteria.
+
 [^Zimmermann2017]:
     > Zimmermann, R.: _A matrix-algebraic algorithm for the Riemannian logarithm on the Stiefel manifold under the canoncial metric.
     > SIAM Journal on Matrix Analysis and Applications 28(2), pp. 322-342, 2017.
@@ -97,13 +106,25 @@ The algorithm is derived in[^Zimmermann2017].
 """
 log(::MetricManifold{ℝ,Stiefel{n,k,ℝ},CanonicalMetric}, ::Any...) where {n,k}
 
+function log(
+    M::MetricManifold{ℝ,Stiefel{n,k,ℝ},CanonicalMetric},
+    p,
+    q;
+    maxiter=1e5,
+    tolerance=1e-9,
+) where {n,k}
+    X = allocate_result(M, log, p, q)
+    log!(M, X, p, q; maxiter=maxiter, tolerance=tolerance)
+    return X
+end
+
 function log!(
     ::MetricManifold{ℝ,Stiefel{n,k,ℝ},CanonicalMetric},
     X,
     p,
     q;
-    tolerance=1e-9,
     maxiter=1e5,
+    tolerance=1e-9,
 ) where {n,k}
     M = p' * q
     QR = qr(q - p * M)
@@ -119,11 +140,13 @@ function log!(
     expnC = exp(-C)
     i = 0
     while (i < maxiter) && (norm(C) > tolerance)
+        println("q\n",q)
         i = i + 1
         LV .= real.(log(V))
         expnC .= exp(-C)
         copyto!(Vpcols, Vpcols * expnC)
     end
+    (i>0) && println(i)
     LV .= real.(LV)
     mul!(X, p, @view(LV[1:k, 1:k]))
     # force the first - Q - to be the reduced form, not the full matrix
diff --git a/test/stiefel.jl b/test/stiefel.jl
index ed82867de9..1b4c8bf9a7 100644
--- a/test/stiefel.jl
+++ b/test/stiefel.jl
@@ -312,6 +312,8 @@ using Manifolds: default_metric_dispatch
         @test isapprox(M3, p, log(M3, p, q), X)
         @test isapprox(M3, p, log(M3, p, r), Y)
         @test inner(M3, p, X, Y) == 0
+        @test inner(M3, p, X, 2*X + 3*Y) == 2*inner(M3, p, X, X)
+        @test norm(M3, p, X) ≈ distance(M3, p, q)
         a = 1.4
         Z = [0.0 2a; -2a 0.0; -a/2 a/2]
         s = exp(M3, p, Z)

From ce20bfbebd2817c9c5597c3ce932d58005bf6b0e Mon Sep 17 00:00:00 2001
From: Ronny <git@ronnybergmann.net>
Date: Thu, 11 Mar 2021 21:00:32 +0100
Subject: [PATCH 21/32] remove some debug print and run formatter.

---
 src/manifolds/StiefelCanonicalMetric.jl | 2 --
 test/stiefel.jl                         | 2 +-
 2 files changed, 1 insertion(+), 3 deletions(-)

diff --git a/src/manifolds/StiefelCanonicalMetric.jl b/src/manifolds/StiefelCanonicalMetric.jl
index 868f8c10bc..aa8e2149f6 100644
--- a/src/manifolds/StiefelCanonicalMetric.jl
+++ b/src/manifolds/StiefelCanonicalMetric.jl
@@ -140,13 +140,11 @@ function log!(
     expnC = exp(-C)
     i = 0
     while (i < maxiter) && (norm(C) > tolerance)
-        println("q\n",q)
         i = i + 1
         LV .= real.(log(V))
         expnC .= exp(-C)
         copyto!(Vpcols, Vpcols * expnC)
     end
-    (i>0) && println(i)
     LV .= real.(LV)
     mul!(X, p, @view(LV[1:k, 1:k]))
     # force the first - Q - to be the reduced form, not the full matrix
diff --git a/test/stiefel.jl b/test/stiefel.jl
index 1b4c8bf9a7..8adce0d0fb 100644
--- a/test/stiefel.jl
+++ b/test/stiefel.jl
@@ -312,7 +312,7 @@ using Manifolds: default_metric_dispatch
         @test isapprox(M3, p, log(M3, p, q), X)
         @test isapprox(M3, p, log(M3, p, r), Y)
         @test inner(M3, p, X, Y) == 0
-        @test inner(M3, p, X, 2*X + 3*Y) == 2*inner(M3, p, X, X)
+        @test inner(M3, p, X, 2 * X + 3 * Y) == 2 * inner(M3, p, X, X)
         @test norm(M3, p, X) ≈ distance(M3, p, q)
         a = 1.4
         Z = [0.0 2a; -2a 0.0; -a/2 a/2]

From 94302ca97a079b52444d9fbbe59214a87825e247 Mon Sep 17 00:00:00 2001
From: kellertuer <git@ronnybergmann.net>
Date: Thu, 18 Mar 2021 09:41:11 +0100
Subject: [PATCH 22/32] improve the test, where an approximation is actually
 needed.

---
 test/stiefel.jl | 13 ++++++++-----
 1 file changed, 8 insertions(+), 5 deletions(-)

diff --git a/test/stiefel.jl b/test/stiefel.jl
index 8adce0d0fb..680027255a 100644
--- a/test/stiefel.jl
+++ b/test/stiefel.jl
@@ -314,10 +314,13 @@ using Manifolds: default_metric_dispatch
         @test inner(M3, p, X, Y) == 0
         @test inner(M3, p, X, 2 * X + 3 * Y) == 2 * inner(M3, p, X, X)
         @test norm(M3, p, X) ≈ distance(M3, p, q)
-        a = 1.4
-        Z = [0.0 2a; -2a 0.0; -a/2 a/2]
-        s = exp(M3, p, Z)
-        Z2 = log(M3, p, s)
-        @test isapprox(M3, p, Z, Z2)
+        # check on a higher dimensional manifold, that the iterations are actually used
+        M4 = MetricManifold(Stiefel(10, 2), CanonicalMetric())
+        p = Matrix{Float64}(I, 10, 2)
+        Random.seed!(42)
+        Z = project(base_manifold(M4), p, randn(size(p)))
+        s = exp(M4, p, Z)
+        Z2 = log(M4, p, s)
+        @test isapprox(M4, p, Z, Z2)
     end
 end

From b13cfb0d3101e001bc3a5381ef0ba1ad3ad1ed8c Mon Sep 17 00:00:00 2001
From: Ronny Bergmann <git@ronnybergmann.net>
Date: Sun, 21 Mar 2021 09:41:03 +0100
Subject: [PATCH 23/32] Update src/manifolds/SymmetricPositiveDefinite.jl

Co-authored-by: Seth Axen <seth.axen@gmail.com>
---
 src/manifolds/SymmetricPositiveDefinite.jl | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/src/manifolds/SymmetricPositiveDefinite.jl b/src/manifolds/SymmetricPositiveDefinite.jl
index 515f00364d..6bff7201bc 100644
--- a/src/manifolds/SymmetricPositiveDefinite.jl
+++ b/src/manifolds/SymmetricPositiveDefinite.jl
@@ -170,4 +170,4 @@ definite matrix `x` on the [`SymmetricPositiveDefinite`](@ref) manifold `M`.
 """
 zero_tangent_vector(::SymmetricPositiveDefinite, ::Any)
 
-zero_tangent_vector!(::SymmetricPositiveDefinite{N}, v, x) where {N} = fill!(v, 0)
+zero_tangent_vector!(::SymmetricPositiveDefinite, v, x) = fill!(v, 0)

From 9b8e3c4ce8cd503eee2e94ea1fa1e3364ae428a7 Mon Sep 17 00:00:00 2001
From: Ronny <git@ronnybergmann.net>
Date: Sun, 21 Mar 2021 09:45:22 +0100
Subject: [PATCH 24/32] address points from code review.

---
 docs/src/manifolds/stiefel.md           | 15 +++++++--------
 src/manifolds/StiefelCanonicalMetric.jl |  8 +++++---
 2 files changed, 12 insertions(+), 11 deletions(-)

diff --git a/docs/src/manifolds/stiefel.md b/docs/src/manifolds/stiefel.md
index a23ef419a5..6227e5c53c 100644
--- a/docs/src/manifolds/stiefel.md
+++ b/docs/src/manifolds/stiefel.md
@@ -19,13 +19,6 @@ Order = [:function]
 ```
 
 ## The canonical metric
-
-```@autodocs
-Modules = [Manifolds]
-Pages = ["manifolds/StiefelCanonicalMetric.jl"]
-Order = [:type]
-```
-
 Any ``X∈T_p\mathcal M``, ``p∈\mathcal M``, can be written as
 
 ```math
@@ -36,9 +29,15 @@ A ∈ ℝ^{p×p} \text{ skew-symmetric},
 B ∈ ℝ^{n×p} \text{ arbitrary.}
 ```
 
-In the Euclidean case, the elements from ``A`` are counted twice (i.e. weighted with a factor of 2).
+In the [`EuclideanMetric`](@ref), the elements from ``A`` are counted twice (i.e. weighted with a factor of 2).
 The canonical metric avoids this.
 
+```@autodocs
+Modules = [Manifolds]
+Pages = ["manifolds/StiefelCanonicalMetric.jl"]
+Order = [:type]
+```
+
 ```@autodocs
 Modules = [Manifolds]
 Pages = ["manifolds/StiefelCanonicalMetric.jl"]
diff --git a/src/manifolds/StiefelCanonicalMetric.jl b/src/manifolds/StiefelCanonicalMetric.jl
index aa8e2149f6..886ae1b05c 100644
--- a/src/manifolds/StiefelCanonicalMetric.jl
+++ b/src/manifolds/StiefelCanonicalMetric.jl
@@ -133,19 +133,21 @@ function log!(
     Vpcols = @view V[:, (k + 1):(2 * k)] #second half of the columns
     S = svd(Vcorner) # preprocessing: Procrustes
     mul!(Vpcols, Vpcols * S.U, S.V')
-    V[:, 1:k] .= [M; QR.R]
+    V[1:k, 1:k] .= M
+    V[(k+1):(2*k), 1:k] .= QR.R
 
     LV = real.(log(V)) # this can be replaced by log_safe.
     C = view(LV, (k + 1):(2 * k), (k + 1):(2 * k))
     expnC = exp(-C)
     i = 0
+    new_Vpcols = Vpcols*expnC # allocate once
     while (i < maxiter) && (norm(C) > tolerance)
         i = i + 1
         LV .= real.(log(V))
         expnC .= exp(-C)
-        copyto!(Vpcols, Vpcols * expnC)
+        mul!(new_Vpcols, Vpcols, expnC)
+        copyto!(Vpcols, new_Vpcols)
     end
-    LV .= real.(LV)
     mul!(X, p, @view(LV[1:k, 1:k]))
     # force the first - Q - to be the reduced form, not the full matrix
     mul!(X, @view(QR.Q[:, 1:k]), @view(LV[(k + 1):(2 * k), 1:k]), true, true)

From 456621e62189e3c74d45ca16d19e5e5f5e735561 Mon Sep 17 00:00:00 2001
From: Ronny <git@ronnybergmann.net>
Date: Sun, 21 Mar 2021 10:08:12 +0100
Subject: [PATCH 25/32] Adress points from code review.

---
 src/manifolds/StiefelCanonicalMetric.jl                 | 8 ++++----
 src/manifolds/SymmetricPositiveDefinite.jl              | 2 +-
 src/manifolds/SymmetricPositiveSemidefiniteFixedRank.jl | 9 ++-------
 3 files changed, 7 insertions(+), 12 deletions(-)

diff --git a/src/manifolds/StiefelCanonicalMetric.jl b/src/manifolds/StiefelCanonicalMetric.jl
index 886ae1b05c..5d6b147f47 100644
--- a/src/manifolds/StiefelCanonicalMetric.jl
+++ b/src/manifolds/StiefelCanonicalMetric.jl
@@ -1,8 +1,8 @@
 @doc raw"""
     CanonicalMetric <: Metric
 
-The Canonical Metric is name used on several manifolds, for example for the [`Stiefel`](@ref)
-manifold, see for example[^EdelmanAriasSmith1998].
+The Canonical Metric refers to a metric for the [`Stiefel`](@ref)
+manifold, see[^EdelmanAriasSmith1998].
 
 [^EdelmanAriasSmith1998]:
     > Edelman, A., Ariar, T. A., Smith, S. T.:
@@ -134,13 +134,13 @@ function log!(
     S = svd(Vcorner) # preprocessing: Procrustes
     mul!(Vpcols, Vpcols * S.U, S.V')
     V[1:k, 1:k] .= M
-    V[(k+1):(2*k), 1:k] .= QR.R
+    V[(k + 1):(2 * k), 1:k] .= QR.R
 
     LV = real.(log(V)) # this can be replaced by log_safe.
     C = view(LV, (k + 1):(2 * k), (k + 1):(2 * k))
     expnC = exp(-C)
     i = 0
-    new_Vpcols = Vpcols*expnC # allocate once
+    new_Vpcols = Vpcols * expnC # allocate once
     while (i < maxiter) && (norm(C) > tolerance)
         i = i + 1
         LV .= real.(log(V))
diff --git a/src/manifolds/SymmetricPositiveDefinite.jl b/src/manifolds/SymmetricPositiveDefinite.jl
index 6bff7201bc..a9daffe529 100644
--- a/src/manifolds/SymmetricPositiveDefinite.jl
+++ b/src/manifolds/SymmetricPositiveDefinite.jl
@@ -170,4 +170,4 @@ definite matrix `x` on the [`SymmetricPositiveDefinite`](@ref) manifold `M`.
 """
 zero_tangent_vector(::SymmetricPositiveDefinite, ::Any)
 
-zero_tangent_vector!(::SymmetricPositiveDefinite, v, x) = fill!(v, 0)
+zero_tangent_vector!(::SymmetricPositiveDefinite, v, ::Any) = fill!(v, 0)
diff --git a/src/manifolds/SymmetricPositiveSemidefiniteFixedRank.jl b/src/manifolds/SymmetricPositiveSemidefiniteFixedRank.jl
index 6b848a09e2..1d8bf5f5a7 100644
--- a/src/manifolds/SymmetricPositiveSemidefiniteFixedRank.jl
+++ b/src/manifolds/SymmetricPositiveSemidefiniteFixedRank.jl
@@ -270,11 +270,6 @@ definite matrix `p` on the [`SymmetricPositiveSemidefiniteFixedRank`](@ref) mani
 """
 zero_tangent_vector(::SymmetricPositiveSemidefiniteFixedRank, ::Any...)
 
-function zero_tangent_vector!(
-    ::SymmetricPositiveSemidefiniteFixedRank{N,K},
-    v,
-    ::Any,
-) where {N,K}
-    fill!(v, 0)
-    return v
+function zero_tangent_vector!(::SymmetricPositiveSemidefiniteFixedRank, v, ::Any)
+    return fill!(v, 0)
 end

From ca261fd26f20a2a1fd5adc3d7c15598c402a30dc Mon Sep 17 00:00:00 2001
From: Seth Axen <seth.axen@gmail.com>
Date: Mon, 22 Mar 2021 00:23:26 -0700
Subject: [PATCH 26/32] Multiply blockwise

---
 src/manifolds/StiefelCanonicalMetric.jl | 5 ++++-
 1 file changed, 4 insertions(+), 1 deletion(-)

diff --git a/src/manifolds/StiefelCanonicalMetric.jl b/src/manifolds/StiefelCanonicalMetric.jl
index 5d6b147f47..aa4bd82ca8 100644
--- a/src/manifolds/StiefelCanonicalMetric.jl
+++ b/src/manifolds/StiefelCanonicalMetric.jl
@@ -60,7 +60,10 @@ function exp!(::MetricManifold{ℝ,Stiefel{n,k,ℝ},CanonicalMetric}, q, p, X) w
     n == k && return mul!(q, p, exp(A))
     QR = qr(X - p * A)
     BC_ext = exp([A -QR.R'; QR.R 0*I])
-    q .= [p Matrix(QR.Q)] * @view(BC_ext[:, 1:k])
+    @views begin
+        mul!(q, p, BC_ext[1:k, 1:k])
+        mul!(q, Matrix(QR.Q), BC_ext[(k + 1):(2 * k), 1:k], true, true)
+    end
     return q
 end
 

From 2b1ebd04c2488d281815f3501a715223dcaf35d4 Mon Sep 17 00:00:00 2001
From: Seth Axen <seth.axen@gmail.com>
Date: Mon, 22 Mar 2021 00:23:38 -0700
Subject: [PATCH 27/32] Use convert to avoid copy

---
 src/manifolds/StiefelCanonicalMetric.jl | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/src/manifolds/StiefelCanonicalMetric.jl b/src/manifolds/StiefelCanonicalMetric.jl
index aa4bd82ca8..b771ce8b90 100644
--- a/src/manifolds/StiefelCanonicalMetric.jl
+++ b/src/manifolds/StiefelCanonicalMetric.jl
@@ -131,7 +131,7 @@ function log!(
 ) where {n,k}
     M = p' * q
     QR = qr(q - p * M)
-    V = Matrix(qr([M; QR.R]).Q * I)
+    V = convert(Matrix, qr([M; QR.R]).Q * I)
     Vcorner = @view V[(k + 1):(2 * k), (k + 1):(2k)] #bottom right corner
     Vpcols = @view V[:, (k + 1):(2 * k)] #second half of the columns
     S = svd(Vcorner) # preprocessing: Procrustes

From 34eba33cac40749889b18c60b9bf0f5954480ef8 Mon Sep 17 00:00:00 2001
From: Seth Axen <seth.axen@gmail.com>
Date: Mon, 22 Mar 2021 00:33:58 -0700
Subject: [PATCH 28/32] Use log_safe!

---
 src/manifolds/StiefelCanonicalMetric.jl | 4 ++--
 1 file changed, 2 insertions(+), 2 deletions(-)

diff --git a/src/manifolds/StiefelCanonicalMetric.jl b/src/manifolds/StiefelCanonicalMetric.jl
index b771ce8b90..b5d3a406e1 100644
--- a/src/manifolds/StiefelCanonicalMetric.jl
+++ b/src/manifolds/StiefelCanonicalMetric.jl
@@ -139,15 +139,15 @@ function log!(
     V[1:k, 1:k] .= M
     V[(k + 1):(2 * k), 1:k] .= QR.R
 
-    LV = real.(log(V)) # this can be replaced by log_safe.
+    LV = log_safe!(allocate(V), V)
     C = view(LV, (k + 1):(2 * k), (k + 1):(2 * k))
     expnC = exp(-C)
     i = 0
     new_Vpcols = Vpcols * expnC # allocate once
     while (i < maxiter) && (norm(C) > tolerance)
         i = i + 1
-        LV .= real.(log(V))
         expnC .= exp(-C)
+        log_safe!(LV, V)
         mul!(new_Vpcols, Vpcols, expnC)
         copyto!(Vpcols, new_Vpcols)
     end

From d7ce92709bbae7178dbe13f6288e7a965d516e08 Mon Sep 17 00:00:00 2001
From: Seth Axen <seth.axen@gmail.com>
Date: Mon, 22 Mar 2021 00:34:08 -0700
Subject: [PATCH 29/32] Avoid unnecessary copy

---
 src/manifolds/StiefelCanonicalMetric.jl | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/src/manifolds/StiefelCanonicalMetric.jl b/src/manifolds/StiefelCanonicalMetric.jl
index b5d3a406e1..ac986edbba 100644
--- a/src/manifolds/StiefelCanonicalMetric.jl
+++ b/src/manifolds/StiefelCanonicalMetric.jl
@@ -146,8 +146,8 @@ function log!(
     new_Vpcols = Vpcols * expnC # allocate once
     while (i < maxiter) && (norm(C) > tolerance)
         i = i + 1
-        expnC .= exp(-C)
         log_safe!(LV, V)
+        expnC = exp(-C)
         mul!(new_Vpcols, Vpcols, expnC)
         copyto!(Vpcols, new_Vpcols)
     end

From 0ceb65b3d2f68ccad50707035b5b0aea13215d50 Mon Sep 17 00:00:00 2001
From: Seth Axen <seth.axen@gmail.com>
Date: Mon, 22 Mar 2021 00:39:30 -0700
Subject: [PATCH 30/32] Multiply blockwise to reduce allocations

---
 src/manifolds/StiefelEuclideanMetric.jl | 8 +++++++-
 1 file changed, 7 insertions(+), 1 deletion(-)

diff --git a/src/manifolds/StiefelEuclideanMetric.jl b/src/manifolds/StiefelEuclideanMetric.jl
index d1d359e7e8..2065da231d 100644
--- a/src/manifolds/StiefelEuclideanMetric.jl
+++ b/src/manifolds/StiefelEuclideanMetric.jl
@@ -30,7 +30,13 @@ exp(::Stiefel, ::Any...)
 
 function exp!(::Stiefel{n,k}, q, p, X) where {n,k}
     A = p' * X
-    return copyto!(q, [p X] * exp([A -X'*X; I A]) * [exp(-A); 0 * I])
+    B = exp([A -X'*X; I A])
+    @views begin
+        r = p * B[1:k, 1:k]
+        mul!(r, X, B[(k + 1):(2 * k), 1:k], true, true)
+    end
+    mul!(q, r, exp(-A))
+    return q
 end
 
 @doc raw"""

From 296f1075a6f1eefeced4bbf2af729224109147c7 Mon Sep 17 00:00:00 2001
From: Seth Axen <seth.axen@gmail.com>
Date: Mon, 22 Mar 2021 00:39:39 -0700
Subject: [PATCH 31/32] Increment version number

---
 Project.toml | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/Project.toml b/Project.toml
index ca6ce76ffc..b2c2ea75eb 100644
--- a/Project.toml
+++ b/Project.toml
@@ -1,7 +1,7 @@
 name = "Manifolds"
 uuid = "1cead3c2-87b3-11e9-0ccd-23c62b72b94e"
 authors = ["Seth Axen <seth.axen@gmail.com>", "Mateusz Baran <mateuszbaran89@gmail.com>", "Ronny Bergmann <manopt@ronnybergmann.net>", "Antoine Levitt <antoine.levitt@gmail.com>"]
-version = "0.4.18"
+version = "0.4.19"
 
 [deps]
 Colors = "5ae59095-9a9b-59fe-a467-6f913c188581"

From 2a2d8e71173fdf3a2236c23c0d2d5c8f8c314748 Mon Sep 17 00:00:00 2001
From: kellertuer <git@ronnybergmann.net>
Date: Mon, 22 Mar 2021 10:35:12 +0100
Subject: [PATCH 32/32] remove windows from testing temporarily.

---
 .github/workflows/ci.yml | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/.github/workflows/ci.yml b/.github/workflows/ci.yml
index 070f5f9ab0..ba20bba9fd 100644
--- a/.github/workflows/ci.yml
+++ b/.github/workflows/ci.yml
@@ -12,7 +12,7 @@ jobs:
     strategy:
       matrix:
         julia-version: ["1.4", "1.5"]
-        os: [ubuntu-latest, macOS-latest, windows-latest]
+        os: [ubuntu-latest, macOS-latest]
     steps:
       - uses: actions/checkout@v2
       - uses: julia-actions/setup-julia@v1