% Source: https://kexue.fm/archives/11864
\documentclass[11pt]{article}
\usepackage[utf8]{inputenc}
\usepackage[T1]{fontenc}
\usepackage{lmodern}
\usepackage{geometry}
\usepackage{hyperref}
\usepackage{graphicx}
\usepackage{amsmath,amssymb}
\usepackage{booktabs,longtable,array}
\usepackage{xcolor}
\usepackage{listings}

\newcommand{\tr}{\mathop{\mathrm{tr}}}
\newcommand{\msign}{\mathop{\mathrm{msign}}}
\newcommand{\skewop}{\mathrm{skew}}

\title{Steepest Descent on Manifolds: 7. Analytical Solution for Stiefel}
\author{Su Jianlin}
\date{August 17, 2026}

\begin{document}

\maketitle

For orthogonal manifolds, our previous conclusion was: in the square matrix case, we can completely write out the analytical solution for the corresponding steepest descent (refer to \href{https://kexue.fm/archives/11215}{\textit{Steepest Descent on Manifolds: 2. Muon + Orthogonal}}); however, for the non-square Stiefel manifold, its steepest descent problem requires solving a system of nonlinear equations, and its analytical solution cannot be written directly (refer to \href{https://kexue.fm/archives/11221}{\textit{Steepest Descent on Manifolds: 3. Muon + Stiefel}}).

But recently, the paper \href{https://papers.cool/arxiv/2608.06218}{\textit{Muon on the Stiefel Manifold Admits an Exact Closed-Form Update}} overturned this assertion, pointing out that an explicit solution can likewise be written for the Stiefel manifold, without needing to solve a system of equations. This article will restate the entire derivation process using my own approach.

\section*{Problem Review}

We will not restate the problem background in too much detail, and directly give the problem to be solved: let $\boldsymbol{W}\in\mathbb{R}^{n\times m}(n \geq m)$ be the current parameters, and $\boldsymbol{G}$ be the gradient of the objective function at $\boldsymbol{W}$. We are looking for the steepest descent direction $\boldsymbol{\Phi}\in\mathbb{R}^{n\times m}$ with step size $\eta$ under the orthogonality constraint, i.e.,
\begin{equation}
\max_{\boldsymbol{\Phi}} \tr(\boldsymbol{G}^{\top}\boldsymbol{\Phi}) \qquad \text{s.t.}\qquad \Vert\boldsymbol{\Phi}\Vert_2 \leq 1,\,\, \boldsymbol{W}^{\top}\boldsymbol{W}=\boldsymbol{I}_m,\,\,(\boldsymbol{W} - \eta \boldsymbol{\Phi})^{\top}(\boldsymbol{W} - \eta \boldsymbol{\Phi})=\boldsymbol{I}_m
\end{equation}
Following the ``first-order approximation is sufficient'' principle, we expand the last constraint and omit the $\eta^2$ term, simplifying it to
\begin{equation}
\max_{\boldsymbol{\Phi}} \tr(\boldsymbol{G}^{\top}\boldsymbol{\Phi}) \qquad \text{s.t.}\qquad \Vert\boldsymbol{\Phi}\Vert_2 \leq 1,\,\, \boldsymbol{W}^{\top}\boldsymbol{W}=\boldsymbol{I}_m,\,\,\boldsymbol{W}^{\top}\boldsymbol{\Phi}+\boldsymbol{\Phi}^{\top}\boldsymbol{W} = \boldsymbol{0}
\end{equation}
where $\Vert\cdot\Vert_2$ is the spectral norm, and $\boldsymbol{I}_m$ is the $m\times m$ identity matrix. Here, changing $\Vert\boldsymbol{\Phi}\Vert_2 \leq 1$ to $\Vert\boldsymbol{\Phi}\Vert_2 = 1$ yields a completely equivalent result (the optimal point of a linear objective function can always be attained on the boundary), but here we retain $\leq$ because $\Vert\boldsymbol{\Phi}\Vert_2 \leq 1$ is a convex set, which provides more room for maneuver during the argument. Our previous solution result was
\begin{equation}
\boldsymbol{\Phi} = \msign(\boldsymbol{G} + \boldsymbol{W}\boldsymbol{X})\label{eq:Phi-ori}
\end{equation}
where $\boldsymbol{X}\in\mathbb{R}^{m\times m}$ is an antisymmetric matrix, and it satisfies the following matrix equation:
\begin{equation}
\boldsymbol{W}^{\top}\msign(\boldsymbol{G} + \boldsymbol{W}\boldsymbol{X})+\msign(\boldsymbol{G} + \boldsymbol{W}\boldsymbol{X})^{\top}\boldsymbol{W} = \boldsymbol{0}\label{eq:X-ori}
\end{equation}
The solution for the case $n=m$ was already obtained in \href{https://kexue.fm/archives/11215}{\textit{Steepest Descent on Manifolds: 2. Muon + Orthogonal}}, but for the standard Stiefel manifold with $n > m$, solving this equation is non-trivial, and even numerical solving is somewhat troublesome. The assertion that ``the Stiefel case has no analytical solution'' arose from this.

\section*{Weakened Version}
\phantomsection\label{sec:weakened}

The key to breaking the deadlock is to remove the equality constraint $\boldsymbol{W}^{\top}\boldsymbol{\Phi}+\boldsymbol{\Phi}^{\top}\boldsymbol{W} = \boldsymbol{0}$ through some explicit representation. To this end, we first introduce a weakened version of the original problem: let $\boldsymbol{\Phi} = \boldsymbol{S}\boldsymbol{W}$, where $\boldsymbol{S}\in\mathbb{R}^{n\times n}$ is any antisymmetric matrix satisfying $\Vert\boldsymbol{S}\Vert_2\leq 1$. Then we consider solving
\begin{equation}
\max_{\boldsymbol{S}} \tr(\boldsymbol{G}^{\top}\boldsymbol{S}\boldsymbol{W}) \qquad \text{s.t.}\qquad \Vert\boldsymbol{S}\Vert_2 \leq 1
\end{equation}
The reason this is a weakening of the original problem is that the $\boldsymbol{\Phi}$ set up in this way obviously satisfies the two constraints of the original proposition:
\begin{align}
1) &\quad \Vert\boldsymbol{\Phi}\Vert_2 = \Vert\boldsymbol{S}\boldsymbol{W}\Vert_2\leq \Vert\boldsymbol{S}\Vert_2 \Vert\boldsymbol{W}\Vert_2\leq 1 \\[5pt]
2) &\quad \boldsymbol{W}^{\top}\boldsymbol{\Phi}+\boldsymbol{\Phi}^{\top}\boldsymbol{W} = \boldsymbol{W}^{\top}\boldsymbol{S}\boldsymbol{W}+\boldsymbol{W}^{\top}\boldsymbol{S}^{\top}\boldsymbol{W} = \boldsymbol{W}^{\top}(\boldsymbol{S} + \boldsymbol{S}^{\top})\boldsymbol{W} = \boldsymbol{0}
\end{align}
So the new problem is solved on a subset of the feasible region of the original problem, and the maximum value obtained will not exceed the maximum value of the original problem. Furthermore, using the trace identity and the antisymmetry of $\boldsymbol{S}$, the objective function can be transformed into
\begin{equation}
\tr(\boldsymbol{G}^{\top}\boldsymbol{S}\boldsymbol{W}) = \tr(\boldsymbol{W}\boldsymbol{G}^{\top}\boldsymbol{S}) = -\tr(\boldsymbol{W}\boldsymbol{G}^{\top}\boldsymbol{S}^{\top}) = -\tr(\boldsymbol{G}\boldsymbol{W}^{\top}\boldsymbol{S}) = \tr([\boldsymbol{W}\boldsymbol{G}^{\top}]_{\skewop}\boldsymbol{S})
\end{equation}
where $[\boldsymbol{X}]_{\skewop} = (\boldsymbol{X} - \boldsymbol{X}^{\top})/2$. So the new problem is equivalent to
\begin{equation}
\max_{\boldsymbol{S}}\tr([\boldsymbol{W}\boldsymbol{G}^{\top}]_{\skewop}\boldsymbol{S}) \qquad \text{s.t.}\qquad \Vert\boldsymbol{S}\Vert_2 \leq 1,\,\,\boldsymbol{S} + \boldsymbol{S}^{\top}=\boldsymbol{0}
\end{equation}
This is already in the form of the standard Muon problem, and from the Muon result we can directly obtain
\begin{equation}
\boldsymbol{S} = \msign([\boldsymbol{G}\boldsymbol{W}^{\top}]_{\skewop})\qquad\Rightarrow\qquad \boldsymbol{\Phi} = \msign([\boldsymbol{G}\boldsymbol{W}^{\top}]_{\skewop})\boldsymbol{W}\label{eq:sol}
\end{equation}
In fact, this is exactly the form of the solution for the square matrix case, and it is also the left rotation in \href{https://kexue.fm/archives/11777}{\textit{Steepest Descent on Manifolds: 6. Muon + Double Rotation}}. One thing worth explaining here is, why do we apply Muon's $\msign$ only after a series of identity transformations? This is mainly to utilize the conclusion that ``the $\msign$ of an antisymmetric matrix is still an antisymmetric matrix'', thereby ensuring that the final result satisfies the antisymmetry requirement of $\boldsymbol{S}$.

\section*{Exact Equivalence}

Next, we will prove that the weakened problem we solved in the previous section actually shares the exact same optimal solution as the original problem!

The proof idea is also very direct. The reason the weakened problem is ``weakened'' is that the current parameterization $\boldsymbol{\Phi} = \boldsymbol{S}\boldsymbol{W}$ is only a subset of the feasible region of the original problem, which might miss some exploration areas and thus fail to reach the optimal solution of the original problem. But if we can prove that for any feasible $\boldsymbol{\Phi}$ in the original problem, we can find an $\boldsymbol{S}$ satisfying the conditions such that $\boldsymbol{\Phi} = \boldsymbol{S}\boldsymbol{W}$, then it shows that their exploration spaces are exactly the same, and thus their optimal solutions are also identical.

The ensuing proof is carried out in three steps. (Note: In hindsight, given appropriate hints, K3 can also independently complete the derivation and proof, refer \href{https://www.kimi.com/share/1a011285-7f82-87e0-8000-0000da55cf07}{here} and \href{https://www.kimi.com/share/1a0110e7-c702-8915-8000-00000fff70c4}{here}.)

\subsection*{Solving the Equation}

The first thing to do is, given $\boldsymbol{\Phi},\boldsymbol{W}$, solve the equation
\begin{equation}
\boldsymbol{\Phi} = \boldsymbol{S}\boldsymbol{W},\qquad \boldsymbol{S}+\boldsymbol{S}^{\top} = \boldsymbol{0}
\end{equation}
Without loss of generality, considering only $n > m$, we can complete $\boldsymbol{W}$ into an $n\times n$ orthogonal matrix $\boldsymbol{U} = [\boldsymbol{W}, \boldsymbol{W}_{\perp}]$. Since $\boldsymbol{S}$ is antisymmetric, $\boldsymbol{U}^{\top}\boldsymbol{S}\boldsymbol{U}$ must also be, so we can write $\boldsymbol{S}$ as
\begin{equation}
\boldsymbol{S} = \boldsymbol{U}\begin{bmatrix}\boldsymbol{A} & -\boldsymbol{B}^{\top} \\
\boldsymbol{B} & \boldsymbol{C}\end{bmatrix}\boldsymbol{U}^{\top} = [\boldsymbol{W}, \boldsymbol{W}_{\perp}]\begin{bmatrix}\boldsymbol{A} & -\boldsymbol{B}^{\top} \\
\boldsymbol{B} & \boldsymbol{C}\end{bmatrix}\begin{bmatrix}\boldsymbol{W}^{\top} \\ \boldsymbol{W}_{\perp}^{\top}\end{bmatrix}\label{eq:S-block}
\end{equation}
where $\boldsymbol{A}\in\mathbb{R}^{m\times m},\boldsymbol{C}\in\mathbb{R}^{(n-m)\times (n-m)}$ are both antisymmetric matrices, and $\boldsymbol{B}\in\mathbb{R}^{(n-m)\times m}$. Then, right-multiplying both sides by $\boldsymbol{W}$ yields
\begin{equation}
\boldsymbol{\Phi} = \boldsymbol{S}\boldsymbol{W} = [\boldsymbol{W}, \boldsymbol{W}_{\perp}]\begin{bmatrix}\boldsymbol{A} & -\boldsymbol{B}^{\top} \\
\boldsymbol{B} & \boldsymbol{C}\end{bmatrix}\begin{bmatrix}\boldsymbol{I}_m \\ \boldsymbol{0}\end{bmatrix} = [\boldsymbol{W}, \boldsymbol{W}_{\perp}]\begin{bmatrix}\boldsymbol{A} \\ \boldsymbol{B}\end{bmatrix}
\end{equation}
Left-multiplying both sides by $\boldsymbol{U}^{\top}$, we can directly read off $\boldsymbol{A} = \boldsymbol{W}^{\top}\boldsymbol{\Phi},\boldsymbol{B} = \boldsymbol{W}_{\perp}^{\top}\boldsymbol{\Phi}$, meaning $\boldsymbol{A},\boldsymbol{B}$ can both be uniquely solved, leaving $\boldsymbol{C}$ as the remaining degree of freedom. Note that $\boldsymbol{A}$ is an antisymmetric matrix, which raises a compatibility condition
\begin{equation}
\boldsymbol{0} = \boldsymbol{A} + \boldsymbol{A}^{\top} = \boldsymbol{W}^{\top}\boldsymbol{\Phi} + \boldsymbol{\Phi}^{\top}\boldsymbol{W}
\end{equation}
which is exactly one of the conditions of the original problem. So far, all results are self-consistent.

\subsection*{Spectral Norm}

The remaining question now is whether we can find an antisymmetric matrix $\boldsymbol{C}$ such that $\Vert \boldsymbol{S}\Vert_2 \leq 1$. The answer is affirmative, but the proof requires a ``\href{https://www.sciencedirect.com/science/article/pii/0022123678900605}{Parrott's Lemma}'' (the original paper cites \href{https://epubs.siam.org/doi/10.1137/0719029}{Davis--Kahan--Weinberger}, but research shows Parrott's lemma is earlier):

\begin{quote}
Given a block matrix $\begin{bmatrix}\boldsymbol{A} & \boldsymbol{C} \\ 
\boldsymbol{B} & \boldsymbol{D}\end{bmatrix}$, fixing $\boldsymbol{A},\boldsymbol{B},\boldsymbol{C}$, we can always find a matrix $\boldsymbol{D}$ such that
\begin{equation}
\left\Vert\begin{bmatrix}\boldsymbol{A} & \boldsymbol{C} \\ 
\boldsymbol{B} & \boldsymbol{D}\end{bmatrix}\right\Vert_2 = \max\left\{\left\Vert\begin{bmatrix}\boldsymbol{A} \\ 
\boldsymbol{B}\end{bmatrix}\right\Vert_2, \left\Vert\begin{bmatrix}\boldsymbol{A} & \boldsymbol{C}\end{bmatrix}\right\Vert_2\right\}
\end{equation}
\end{quote}

Returning to our problem. To make $\Vert \boldsymbol{S}\Vert_2 \leq 1$, we only need $\left\Vert\begin{bmatrix}\boldsymbol{A} & -\boldsymbol{B}^{\top} \\
\boldsymbol{B} & \boldsymbol{C}\end{bmatrix}\right\Vert_2\leq 1$. For this block matrix, we have
\begin{equation}
\left\Vert\begin{bmatrix}\boldsymbol{A} \\
\boldsymbol{B}\end{bmatrix}\right\Vert_2 = \left\Vert\begin{bmatrix}\boldsymbol{W}^{\top}\boldsymbol{\Phi} \\
\boldsymbol{W}_{\perp}^{\top}\boldsymbol{\Phi}\end{bmatrix}\right\Vert_2 = \Vert\boldsymbol{U}^{\top}\boldsymbol{\Phi}\Vert_2 = \Vert\boldsymbol{\Phi}\Vert_2 \leq 1
\end{equation}
As for $[\boldsymbol{A},-\boldsymbol{B}^{\top}] = [-\boldsymbol{A}^{\top},-\boldsymbol{B}^{\top}]$, which is the negative transpose of $\begin{bmatrix}\boldsymbol{A} \\
\boldsymbol{B}\end{bmatrix}$, the spectral norm remains unchanged (and likewise does not exceed 1). Therefore, according to Parrott's lemma, we can find a matrix $\boldsymbol{C}$ such that $\Vert \boldsymbol{S}\Vert_2 = \left\Vert\begin{bmatrix}\boldsymbol{A} & -\boldsymbol{B}^{\top} \\
\boldsymbol{B} & \boldsymbol{C}\end{bmatrix}\right\Vert_2 \leq 1$.

\subsection*{Antisymmetry}

However, things are not over yet. We also have an antisymmetry requirement for $\boldsymbol{C}$, which Parrott's lemma does not guarantee.

Fortunately, this patch is not difficult. Let the matrix given by Parrott's lemma be $\boldsymbol{C}_0$, we directly take $\boldsymbol{C} = [\boldsymbol{C}_0]_{\skewop}$. Since the remaining parts inherently satisfy the antisymmetry requirement, we have
\begin{equation}
\begin{bmatrix}\boldsymbol{A} & -\boldsymbol{B}^{\top} \\
\boldsymbol{B} & [\boldsymbol{C}_0]_{\skewop}\end{bmatrix} = \frac{1}{2}\left(\begin{bmatrix}\boldsymbol{A} & -\boldsymbol{B}^{\top} \\
\boldsymbol{B} & \boldsymbol{C}_0\end{bmatrix} - \begin{bmatrix}\boldsymbol{A} & -\boldsymbol{B}^{\top} \\
\boldsymbol{B} & \boldsymbol{C}_0\end{bmatrix}^{\top}\right) 
\end{equation}
The key point here is that antisymmetrization does not increase the spectral norm: by the triangle inequality, we get
\begin{equation}
\left\Vert\begin{bmatrix}\boldsymbol{A} & -\boldsymbol{B}^{\top} \\
\boldsymbol{B} & [\boldsymbol{C}_0]_{\skewop}\end{bmatrix}\right\Vert_2 \leq \frac{1}{2}\left\Vert\begin{bmatrix}\boldsymbol{A} & -\boldsymbol{B}^{\top} \\
\boldsymbol{B} & \boldsymbol{C}_0\end{bmatrix}\right\Vert_2 + \frac{1}{2}\left\Vert \begin{bmatrix}\boldsymbol{A} & -\boldsymbol{B}^{\top} \\
\boldsymbol{B} & \boldsymbol{C}_0\end{bmatrix}^{\top}\right\Vert_2 \leq \frac{1}{2} + \frac{1}{2} = 1 
\end{equation}
So $[\boldsymbol{C}_0]_{\skewop}$ is the antisymmetric matrix we are ultimately looking for.

At this point, we have proved the result we initially hoped to prove: for every $\boldsymbol{\Phi}$ in the original problem, we can find an $n\times n$, antisymmetric matrix $\boldsymbol{S}$ with a spectral norm not exceeding 1, such that $\boldsymbol{\Phi}=\boldsymbol{S}\boldsymbol{W}$. Thus, the exploration space of the weakened version is identical to the exploration space of the original problem. In other words, Equation \eqref{eq:sol} is the exact solution to the original problem!

\section*{Efficient Computation}

From a theoretical perspective, the problem has been perfectly resolved; but from a practical perspective, there is still an efficiency issue worth discussing. Muon requires performing $\msign$ on an $n\times m$ matrix, whereas the exact solution \eqref{eq:sol} on the Stiefel manifold requires performing $\msign$ on an $n\times n$ matrix. If computed directly, it is obviously much more expensive than Muon when $n \gg m$, so it is necessary to explore schemes that save computation.

Actually, it is not hard to imagine. First, we rewrite the core operation into a low-rank form:
\begin{equation}
\msign([\boldsymbol{G}\boldsymbol{W}^{\top}]_{\skewop}) = \msign(\boldsymbol{G}\boldsymbol{W}^{\top} - \boldsymbol{W}\boldsymbol{G}^{\top}) = \msign\bigg(\begin{bmatrix}\boldsymbol{G} & \boldsymbol{W}\end{bmatrix} \underbrace{\begin{bmatrix}\boldsymbol{0} & \boldsymbol{I}_m \\ -\boldsymbol{I}_m & \boldsymbol{0}\end{bmatrix}}_{\boldsymbol{J}} \begin{bmatrix}\boldsymbol{G}^{\top} \\ \boldsymbol{W}^{\top}\end{bmatrix}\bigg)
\end{equation}
It can be seen that the rank of $[\boldsymbol{G}\boldsymbol{W}^{\top}]_{\skewop}$ is at most $2m$. Assuming $n \geq 2m$, we first perform a QR decomposition on $[\boldsymbol{G},\boldsymbol{W}]$ to get $\boldsymbol{Q}\boldsymbol{R}$, where $\boldsymbol{Q}\in\mathbb{R}^{n\times 2m}$ satisfies $\boldsymbol{Q}^{\top}\boldsymbol{Q}=\boldsymbol{I}_{2m}$, and $\boldsymbol{R}\in\mathbb{R}^{2m\times 2m}$ is an upper triangular matrix. Then the part to be calculated becomes $\msign(\boldsymbol{Q}\boldsymbol{R}\boldsymbol{J}\boldsymbol{R}^{\top} \boldsymbol{Q}^{\top})$. Note that the $\msign$ operation has covariance with respect to orthogonal matrices, i.e.,
\begin{equation}
\msign(\boldsymbol{Q}\boldsymbol{R}\boldsymbol{J}\boldsymbol{R}^{\top} \boldsymbol{Q}^{\top}) = \boldsymbol{Q}\msign(\boldsymbol{R}\boldsymbol{J}\boldsymbol{R}^{\top} )\boldsymbol{Q}^{\top}
\end{equation}
So in reality, we only need to perform $\msign$ once on the $2m\times 2m$ sized matrix $\boldsymbol{R}\boldsymbol{J}\boldsymbol{R}^{\top}$, which provides a significant speedup for scenarios where $n\gg 2m$. However, if $n$ and $m$ do not differ by more than an order of magnitude, it is more convenient to calculate using the original formula, since QR decomposition also incurs a cost.

In addition, how to switch to a momentum scenario is also worth discussing. We have two conceivable schemes: the first is $\boldsymbol{M} = \mathrm{EMA}(\boldsymbol{G})$, in which case we only need to replace $\boldsymbol{G}$ with $\boldsymbol{M}$, leaving the rest unchanged; the other is $\boldsymbol{M} = \mathrm{EMA}(\boldsymbol{G}\boldsymbol{W}^{\top})$, in which case $[\boldsymbol{G}\boldsymbol{W}^{\top}]_{\skewop}$ is replaced by $[\boldsymbol{M}]_{\skewop}$. This form can retain more information, at the cost of occupying more space. Moreover, after a long-term moving average, even if $n\gg 2m$, $\boldsymbol{M}$ typically will not be low-rank, and thus cannot be accelerated.

Of course, under mainstream model architectures, scenarios where $n\gg 2m$ are not common (only in the case of Per-Head Muon, the matrix for each Head barely qualifies), so the practical significance of this acceleration trick is limited, and we can just briefly understand it.

\section*{Open Problem}

Here the author raises an open question: if we remove the orthogonality constraint and only consider steepest descent on the tangent manifold, what would its optimal solution be? That is, considering
\begin{equation}
\max_{\boldsymbol{\Phi}} \tr(\boldsymbol{G}^{\top}\boldsymbol{\Phi}) \qquad \text{s.t.}\qquad \Vert\boldsymbol{\Phi}\Vert_2 \leq 1,\,\,\boldsymbol{W}^{\top}\boldsymbol{\Phi}+\boldsymbol{\Phi}^{\top}\boldsymbol{W} = \boldsymbol{0}
\end{equation}
This problem can be viewed as a generalization of ``Stiefel-version Muon'' to general matrices; it does not depend on the orthogonality of $\boldsymbol{W}$ itself, and thus is expected to be applicable to more general scenarios. The result the author can currently obtain is that Equation \eqref{eq:Phi-ori} and Equation \eqref{eq:X-ori} are still applicable, meaning the optimal solution has the following form:
\begin{equation}
\boldsymbol{\Phi} = \msign(\boldsymbol{G} + \boldsymbol{W}\boldsymbol{X})
\end{equation}
where $\boldsymbol{X}\in\mathbb{R}^{m\times m}$ is an antisymmetric matrix satisfying the following matrix equation:
\begin{equation}
\boldsymbol{W}^{\top}\msign(\boldsymbol{G} + \boldsymbol{W}\boldsymbol{X})+\msign(\boldsymbol{G} + \boldsymbol{W}\boldsymbol{X})^{\top}\boldsymbol{W} = \boldsymbol{0}
\end{equation}
As for how to solve it when $\boldsymbol{W}^{\top}\boldsymbol{W}\neq\boldsymbol{I}_m$, it remains unknown for the time being (assuming $n=m$ seems unable to simplify it either). Furthermore, the construction of the \hyperref[sec:weakened]{weakened problem} also holds for it, requiring only a slight modification:
\begin{equation}
\boldsymbol{\Phi}_{\text{weak}} = \msign([\boldsymbol{G}\boldsymbol{W}^{\top}]_{\skewop})\boldsymbol{W}/\Vert\boldsymbol{W}\Vert_2
\end{equation}
However, because the earlier equivalence proof strongly depends on $\boldsymbol{W}^{\top}\boldsymbol{W}=\boldsymbol{I}_m$, it appears that this weak solution is genuinely just a weak solution, rather than the optimal solution. Everyone is welcome to provide more progress~

\section*{Article Summary}

In this article, we successfully found an analytical solution for the manifold steepest descent of ``Muon + Stiefel'', bringing this problem to a satisfactory conclusion. From now on, regardless of whether it is a square matrix or not, the Muon optimizer under orthogonal constraints can be calculated analytically, and solving it itself will no longer be a difficulty.

\vspace{30pt}
\noindent\textcolor[HTML]{FF8800}{\textit{\textbf{When reprinting, please include the address of this article:} \href{https://kexue.fm/archives/11864}{https://kexue.fm/archives/11864}}}

\vspace{15pt}
\noindent\textcolor[HTML]{FF8800}{\textit{\textbf{For more detailed reprinting matters, please refer to:} \href{https://kexue.fm/archives/6508}{\textit{Scientific Space FAQ}}}}

\end{document}
