\section{Introduction}
	In this paper, we consider the following nonsmooth optimization problem
	\begin{equation}
		\label{Prob_Ori}
		\tag{OCP}
		\begin{aligned}
			\min_{X \in \bb{R}^{n\times p}} \quad &f(X)\\
			\text{s. t.} \quad & X\tp X = I_p,
		\end{aligned}
	\end{equation}
	where the objective function $f: \bb{R}^{n\times p} \mapsto \bb{R}$ is locally Lipschitz continuous and possibly nonsmooth, while  $I_p$ denotes the $p\times p$ identity matrix.
	We denote the feasible region of \ref{Prob_Ori} as $ \ca{S}_{n,p} := \{X \in \bb{R}^{n\times p}:  X\tp  X = I_p \}$, which is also referred as the Stiefel manifold throughout this paper.

	Our interest in \ref{Prob_Ori} with nonsmooth objective functions comes from its various applications in different engineering fields, including sparse principal component analysis \cite{chen2018proximal,xiao2020l21}, robust subspace recovery \cite{lerman2018overview,wang2019globally},  packing problems \cite{absil2017collection}, and  training orthogonally constrained  neural networks  \cite{arjovsky2016unitary,bansal2018can}, where the objective function is only locally Lipschitz continuous and usually not weakly convex (i.e. $f(X) + \frac{\tau}{2} \norm{X}\fs$ is convex for some constant $\tau \geq 0$ \cite{davis2019stochastic}). As an illustration, let us present some motivating applications arisen from training deep neural networks.

	Training deep neural networks is usually thought to be challenging both theoretically and practically, for which the vanishing/exploding gradients is one of the most important reasons \cite{glorot2010understanding}.  To address such issue, several recent works focus on imposing orthogonality constraints to the weights of the layers in these deep neural networks \cite{arjovsky2016unitary,bansal2018can,lezcano2019trivializations}. As orthogonality implies energy preservation properties  \cite{zhou2006special}, these existing works demonstrate that the orthogonal constraints can stabilize the distribution of activations over layers within convolutional neural networks and make their optimization more efficient. Moreover, some existing works \cite{lezcano2019cheap,lezcano2019trivializations,wang2020orthogonal} observe encouraging improvements in the accuracy and robustness of the networks with orthogonal constraints.
	However, when these neural networks are built from nonsmooth activation functions, their loss functions are usually not weakly convex.   For example, the loss function for fully connected two-layer neural network with mean squared error (MSE) can be expressed as
	\begin{equation*}
		f_{\text{mse}(W_1, W_2, b_1, b_2)\,=\,}\frac{1}{N}\sum_{i = 1}^{N} ~ \norm{ \tilde{\sigma}_2(W_2\tilde{\sigma}_1(W_1x_i + b_1) + b_2) - y_i }_2^2,
	\end{equation*}
	where $\{W_1, W_2\}$ denote the weight matrices, $\{b_1, b_2\}$ denote the bias vectors, $\{\tilde{\sigma}_1, \tilde{\sigma}_2\}$ refer to the activation functions, and $\{(x_i, y_i)\mid i = 1,...,N\}$ is the training data set.
	When $\tilde{\sigma}_i$ ($i=1,2$) are chosen as  nonsmooth activation functions, such as  rectified linear unit (ReLU) or leaky ReLU \cite{maas2013rectifier}, the loss function is clearly locally Lipschitz continuous but not weakly convex with respect to $(W_1, W_2, b_1, b_2)$.
	Therefore, we only assume that the objective function $f$ in \ref{Prob_Ori} to be locally Lipschitz continuous throughout this paper.

	\subsection{Existing works}
	Minimizing smooth objective functions over the Stiefel manifold has been extensively studied in the past several years.
	Interested readers can refer to the books \cite{Absil2009optimization,boumal2020introduction}, the recent survey paper \cite{hu2020brief} and the references therein for more details.
	Compared with smooth optimization over the Stiefel manifold, research on nonsmooth cases are limited \cite{li2019weakly}.
	In the following, we briefly mention some existing state-of-the-art approaches for nonsmooth optimization over the Stiefel manifold.

	\paragraph{Riemannian subgradient methods}

	Some existing Riemannian subgradient methods are specifically designed for minimizing geodsically convex objectives over a Riemannian manifold \cite{ferreira1998subgradient,ferreira2002proximal,zhang2016first}. However, the geodesic convexity over a compact Riemannian manifold, such as the Stiefel manifold discussed in this paper, only holds for constant functions. Therefore, these methods are invalid for any nontrivial case of \ref{Prob_Ori}.
	Very recently, \cite{li2019weakly} studies a class of Riemannian subgradient methods to minimize weakly convex functions over the Stiefel manifold. By extending the subgradient inequality \cite{davis2019stochastic} from $\bb{R}^{n\times p}$ to the tangent spaces of the Stiefel manifold, the theoretical analysis of these Riemannian subgradient methods can be implemented by the same methodologies as existing works for unconstrained optimization \cite{davis2018subgradient,davis2019stochastic,li2019incremental}. However, the analysis in \cite{li2019weakly} relies on the weak convexity of the objective function.  Therefore, when we apply subgradient methods to solve \ref{Prob_Ori}, the proposed frameworks in \cite{li2019weakly} are not capable of analyzing the related theoretical properties.

	\paragraph{Riemannian gradient sampling methods}

	Gradient sampling methods \cite{burke2005robust,burke2020gradient}
	are originally proposed for unconstrained nonsmooth nonconvex optimization. For nonsmooth optimization problems over a Riemannian manifold, several recent works \cite{hosseini2018line,hosseini2017riemannian} proposed Riemannian gradient sampling methods based on their unconstrained origins. When applied to the minimization over the Stiefel manifold, those Riemannian gradient sampling methods first take a few sampling points $\{ X_{k,j}\mid j = 1,...,N_J\} \subset \ca{S}_{n,p}$ in a neighborhood of the current iterate $X_k$. Then they aim to find a descent direction in the convex hull of the Riemannian gradient at points $\{ X_{k,j} \mid j = 1,...,N_J\}$. However, as mentioned in \cite{hosseini2017riemannian,hosseini2018line}, the number of sample points $N_J$ should be much larger than the dimension of the manifold. Therefore, for optimization over the Stiefel manifold, the  number of sampling points should be much larger than $np-p(p+1)/2$ as described in \cite{li2019weakly,burke2020gradient}. As emphasized in \cite{li2019weakly}, it is usually expensive to generate a descent direction in existing Riemannian gradient sampling methods, especially in high dimensional cases.

	\paragraph{Proximal gradient methods}

	For problems with prox-friendly objective functions, for instance, the summation of a
	smooth function and a nonsmooth regularizer whose proximal mapping is easy-to compute,
	the Riemannian proximal gradient methods are developed by extending the proximal gradient methods from Euclidean space to Riemannian manifolds. Several existing works \cite{ferreira2002proximal,de2016new} develop the Riemannian proximal gradient methods by computing the proximal mappings over the Riemannian manifold. Although their global convergence properties could be established by following the same techniques as their unconstrained counterparts, computing the proximal mappings in these approaches is as difficult as the original problem, hence they are generally inefficient in practice \cite{li2019weakly}. Recently, several Riemannian proximal gradient methods \cite{chen2018proximal,huang2019riemannian,zhou2021semi} are proposed by computing the proximal mapping in the tangent space rather than over the manifold. Therefore, computing the proximal mappings in each iteration is equivalent to solving a linearly constrained strongly convex optimization problem, for which various existing efficient solvers such as semi-smooth Newton methods \cite{qi1993nonsmooth} and Arrow-Hurwicz methods \cite{chambolle2011a} can be applied.

	Apart from these Riemannian proximal gradient methods, a variety of existing approaches are developed by regarding \ref{Prob_Ori} as a constrained optimization problem in $\bb{R}^{n\times p}$. Some of these approaches are built upon splitting and alternating techniques, such as the splitting method for orthogonality constrained problem (SOC) \cite{lai2014splitting}, ALM-based variable splitting framework \cite{rosman2014augmented-lagrangian}, proximal alternating minimization approach based on the augmented Lagrangian method (PAMAL) \cite{chen2016augmented}, etc. By splitting the objective function and constraints apart properly, the subproblems in these approaches are usually easy to compute. However, these approaches usually lack convergence guarantees and their performance is very sensitive to the parameters, which are difficult to tune in practice \cite{chen2018proximal,li2019weakly,huang2019riemannian}.

	Recently, a novel class of efficient approaches for Stiefel manifold optimization are developed based on  exact penalty models. Inspired by the exact penalty model (PenC) for smooth optimization over the Stiefel manifold \cite{xiao2020class,hu2020anefficiency}, \cite{xiao2020l21} extends PenC to $\ell_{2,1}$-norm regularized cases and proposes a proximal gradient method called PenCPG.
	In PenCPG, the proximal subproblem has a closed-form solution, which leads
	to its numerical superiority over existing Riemannian proximal gradient approaches.

	Furthermore, for generalized composite objective functions, \cite{xiao2021penalty}
	proposed a penalty-free infeasible approach named sequential linearized proximal gradient method (SLPG).
	In each iteration, SLPG alternatively takes the tangential and the normal steps, both of which do not involve any
	penalty parameter or orthonormalization procedure. Consequently, \SLPGs enjoys high scalability and avoids the numerical inefficiency from inappropriately selected penalty parameters.

	However, the efficiencies of all the aforementioned proximal gradient methods heavily rely on having
	prox-friendly objective functions.
	For general nonsmooth objective functions, computing its corresponding proximal mappings can be extremely expensive \cite{li2019weakly}. As a result, we need new approaches for general nonsmooth optimization over the Stiefel manifold.

	\subsection{Motivation}

	Our motivation comes from the penalty function approaches for Riemannian optimization. For smooth optimization problems on the Stiefel manifold, \cite{xiao2021solving} proposes an exact penalty  function (ExPen) that takes the form as follows:
	\begin{equation}
		\tag{ExPen}
		\label{ExPen}
		f\left( X\Apen{X}\right) + \frac{\beta}{4} \norm{X\tp X - I_p}\fs.
	\end{equation}
	They show that \ref{ExPen} is an exact penalty function for smooth optimization over the Stiefel manifold. Different from the existing Fletcher's penalty function \cite{fletcher1970class}, the objective function in \ref{ExPen} does not involve $\nabla f$. Therefore, \ref{ExPen} has easy-to-compute derivatives and enables direct implementation of various existing unconstrained approaches for solving smooth optimization problems over the Stiefel manifold.

	Very recently,  \cite{xiao2022constraint} proposes constraint dissolving approaches for minimizing smooth functions over a closed Riemannian manifold. In their proposed approaches, solving a Riemannian optimization problem is transferred into  the unconstrained minimization of a corresponding constraint dissolving function. According to \cite{xiao2022constraint}, the  constraint dissolving function for \ref{Prob_Ori} can be expressed as
	\begin{equation}
		\tag{CDF}
		\label{CDF}
		f\left( \A(X) \right) + \frac{\beta}{4} \norm{X\tp X - I_p}\fs,
	\end{equation}
	where the constraint dissolving mapping $\A: \bb{R}^{n\times p} \to \bb{R}^{n\times p}$ is locally Lipschitz smooth and  satisfies the following conditions.

	\begin{cond}\label{cond:1}
		\begin{enumerate}
			\item $\A(X) = X$ holds for any $X \in \ca{S}_{n,p}$.
			\item The Jacobian of $\A(X)\tp \A(X) - I_p$ equals to $0$ for any $X \in \ca{S}_{n,p}$.
		\end{enumerate}
	\end{cond}

	Condition \ref{cond:1}
	grants great flexibility in designing the constraint dissolving functions for \ref{Prob_Ori}.
	However, for nonsmooth optimization problems over the Stiefel manifold, the chain rule usually fails when computing $\partial (f\circ \A)$. Therefore, the Clarke subdifferential of \ref{CDF} is usually not easy to achieve  from $\partial f$ and the Jacobian of $\A$.
	Meanwhile, a sufficient condition to guarantee the validity of chain rule is that $\A$ is a {\it homeomorphism} (i.e. $\A^{-1}$ is well-defined and continuous over $\bb{R}^{n\times p}$) \cite[Theorem 2.3.10]{clarke1990optimization}.

	\subsection{Contribution}
	Taking both Condition \ref{cond:1} and the desirable homeomorphism property of $\A$ into account,
	we construct the following nonsmooth constraint dissolving function named \ref{NEPen}
	for nonsmooth optimization over the Stiefel manifold,
	\begin{equation}
		\label{NEPen}
		\tag{NCDF}
		h(X):= f\left( \A(X) \right) + \frac{\beta}{4} \norm{X\tp X - I_p}\fs,
	\end{equation}
	where $\A: \bb{R}^{n\times p} \to \bb{R}^{n\times p}$ is defined as
	\begin{equation}
		\A(X) := \frac{1}{8}X \left( {15} I_p - {10}X\tp X + 3 (X\tp X)^2  \right).
	\end{equation}
	We prove that $\A$ is a homeomorphism and $\partial (f\circ \A)$ follows the chain rule \cite[Theorem 2.3.10]{clarke1990optimization},  and  $\partial h(X)$ has an explicit expression
	which is easy to achieve from the subdifferential of $f$ . Based on the explicit expression of $\partial h(X)$, we prove that  \ref{Prob_Ori} and \ref{NEPen} share the same local minimizers and stationary points under mild assumptions. These properties characterize the equivalence between \ref{Prob_Ori} and \ref{NEPen} and further illustrate that various unconstrained approaches can be directly implemented to solve \ref{Prob_Ori} through \ref{NEPen}.

	Moreover, we present several examples to show that \ref{Prob_Ori} can be solved by directly applying subgradient methods to minimize \ref{NEPen}. In these examples, we propose a class of subgradient methods and show that their convergence properties can be established directly from existing rich theoretical results in unconstrained optimization. Preliminary numerical experiments on training neural networks with normalized weights further illustrate that \ref{NEPen} can be minimized by direct and efficient implementations of various existing advanced unconstrained optimization solvers.  These illustrative examples highlight the great potential of \ref{NEPen}.

	\subsection{Notations}
	Throughout this paper, the Euclidean inner product of two matrices $X, Y\in \bb{R}^{n\times p}$ is defined as $ \inner{X, Y}=\tr(X\tp Y)$,
	where $\tr(A)$ is the trace of a matrix $A\in \bb{R}^{n\times p}$.
	We use $\norm{\cdot}_2$ and $\norm{\cdot}\ff$ to represent the $2$-norm and the Frobenius norm, respectively.
	The notations $\mathrm{diag}(A)$ and $\Diag(x)$
	stand for the vector formed by the diagonal entries of the matrix $A$,
	and the diagonal matrix with the entries of $x\in\bb{R}^n$ as its diagonal, respectively.
	We denote the smallest eigenvalue of a square matrix $A$ by $\lambda_{\mathrm{min}}(A)$, and set $$\Phi(M):= \frac{1}{2}(M+M\tp).$$
	Moreover,  we set the metric on the Stiefel manifold as the metric inherited from
	the inner product in $\bb{R}^{n\times p}$. Then we denote
	$\ca{T}_X$ as the tangent space of the
	Stiefel manifold at $X$, which can be expressed as
	\begin{equation*}
		\ca{T}_X := \{D \in \bb{R}^{n\times p}: \Phi(D\tp X) = 0 \},
	\end{equation*}
	while $\ca{N}_X$ denotes the normal space of the Stiefel manifold at $X$,
	\begin{equation*}
		\ca{N}_X := \{D \in \bb{R}^{n\times p}: D = X\Lambda, \Lambda = \Lambda\tp \}.
	\end{equation*}
	Additionally, $\ProjS(X) = UV\tp$ denotes the orthogonal projection to Stiefel manifold, where $X = U\Sigma V\tp$ is the economical SVD of $X$
	with $U\in \ca{S}_{n,p}$, $V\in \ca{S}_{p,p}$ and $\Sigma$
	is a $p\times p$ diagonal matrix with the singular values of
	$X$ on its diagonal.  Moreover, the ball centered at $X$ with radius $\delta$ is defined as $\ca{B}(X, \delta):=  \{ Y \in \bb{R}^{n\times p}: \norm{Y - X}\ff \leq \delta  \}$.
	Furthermore, $\JA(X)[D]$ refers to the linear transform of $D$ by the linear mapping $\JA(X)$, i.e., $\JA(X)[D] = \lim_{t \to 0} \frac{1}{t} (\A(X+tD) - \A(X))$. Additionally, for any subset $\ca{F} \subseteq \bb{R}^{n\times p}$, we denote
	\begin{equation*}
		\JA(X)[\ca{F}] := \{ \JA(X)[D]: D \in \ca{F} \}.
	\end{equation*}
	Finally, we denote $g(X) := f(\A(X))$, and
	$$
	 \Omega_r = \Big\{ X \in \bb{R}^{n\times p} \,:\, \norm{X^\top X - I_p}\ff \leq r
	 \Big\}.
	$$

	\subsection{Organization}
	The rest of this paper is organized as follows. Section 2 presents several preliminary notations and definitions.  We analyze the relationships between \ref{Prob_Ori} and \ref{NEPen} in Section 3. We show how to implement subgradient methods for \ref{Prob_Ori} by \ref{NEPen} and prove its theoretical convergence directly from existing works in Section 4. In Section 5, we present preliminary numerical experiments to illustrate that \ref{NEPen} enables efficient and direct implementation of various existing unconstrained solvers to solve \ref{Prob_Ori}. We draw a brief conclusion in the last section.

	\section{Preliminaries}
	\subsection{Definitions}

	\begin{defin}

		Given $X \in \bb{R}^{n\times p}$, the generalized directional derivative of $f$ at $X$  in the direction $D \in \bb{R}^{n\times p}$, denoted by $f^\circ(X, D)$, is defined as
		\begin{equation*}
			f^\circ(X, D) = \mathop{\lim\sup}_{Y\to X, t \downarrow  0}~ \frac{f(Y +tD) - f(Y)}{t}.
		\end{equation*}
		Then the generalized gradient or the Clarke subdifferential of $f$ at $X \in \bb{R}^{n\times p}$, denoted by $\partial f(X)$, is defined as
		\begin{equation*}
			\partial f(X) := \{ W \in \bb{R}^{n\times p} : \inner{W, D} \leq f^\circ (X, D), \text{ for all } D \in \bb{R}^{n\times p} \}.
		\end{equation*}
	\end{defin}

	\begin{defin}

		We say that $f$ is regular at $X \in \bb{R}^{n\times p}$ if for every direction $D$, the one-sided directional derivative
		\begin{equation*}
			f^\star(X,D) = \lim_{t \downarrow 0} \frac{f(X + tD) - f(X)}{t}
		\end{equation*}
		exists and $f^\star(X, D) = f^\circ(X, D)$.
	\end{defin}

	Next we follow the definition in \cite{yang2014optimality,hosseini2018line} to present the definition for Riemannian subdifferential on the Stiefel manifold.
	\begin{defin}
		Given $X \in \ca{S}_{n,p}$, the generalized Riemannian directional derivative of $f$ at $X$ in the direction $D \in \ca{T}_X$, denoted by $f^\diamond(X, D)$, is given by
		\begin{equation*}
			f^\diamond(X, D) = \mathop{\lim\sup}_{ \ca{S}_{n,p} \ni Y\to X, t \downarrow  0}~ \frac{f\left( \ProjS(Y + tD) \right) - f(Y)}{t}.
		\end{equation*}
		Then the Riemannian subdifferential, denoted as $\partialRie f(X)$, is defined as
		\begin{equation*}
			\partialRie f(X) := \{ W \in \ca{T}_X : \inner{W, D} \leq f^\diamond (X, D), ~\text{for all } D \in  \ca{T}_X \}.
		\end{equation*}
	\end{defin}

	\begin{defin}
		\label{Defin_Projected_subdifferential}
		Given any $X\in \ca{S}_{n,p}$, the projected subdifferential of $f$ at $X \in \ca{S}_{n,p}$ on the Stiefel manifold is defined as
		\begin{equation*}
			\partialProj f(X) = \{ W - X\Phi(X\tp W): W \in \partial f(X) \}.
		\end{equation*}
	\end{defin}

	Moreover, we define the first-order stationary points of problem \ref{Prob_Ori} as follows.
	\begin{defin}[{\cite[Theorem 6.11]{clarke1990optimization}}]
		\label{Defin_FOSP}
		Given $X \in \ca{S}_{n,p}$, we say $X$ is a first-order stationary point of \ref{Prob_Ori} if and only if
		\begin{equation*}
			0 \in \partialProj f(X).
		\end{equation*}
	\end{defin}

	\begin{prop}[{\cite[Theorem 5.1]{yang2014optimality}}]
		\label{Prop_relationship_partialR_partialP}
		For any given $X \in \ca{S}_{n,p}$, it holds that
		\begin{equation*}
			\partialRie f(X) \subseteq \partialProj f(X).
		\end{equation*}
		Furthermore,  the equality holds when $f$ is Clarke regular.
	\end{prop}

	\begin{rmk}
		It is worth mentioning that Proposition \ref{Prop_relationship_partialR_partialP} implies that any local minimizer $\tX$ of \ref{Prob_Ori} satisfies $0 \in \partialProj f(X)$.  Moreover,  Definition \ref{Defin_FOSP} coincides with the widely used optimality conditions \cite{yang2014optimality} for regular objective functions. Additionally, since $\partial f(X)$ is usually easy to compute, the projected subdifferential $\partialProj f(X)$ is usually much easier to compute than $\partialRie f(X)$.
	\end{rmk}

	\begin{defin}[\cite{clarke1990optimization,rockafellar2009variational}]
		For any locally Lipschitz continuous function $\psi: \bb{R}^{n\times p} \to \bb{R}$, we say $X \in \bb{R}^{n\times p}$ is a first-order stationary point of function $\psi$ if and only if
		\begin{equation*}
			0 \in \partial \psi(X).
		\end{equation*}
	\end{defin}

	The following theorem is a direct corollary from \cite[Theorem 2.3.10]{clarke1990optimization} that characterizes $\partial g(X)$, where $g := f\circ \A.$
	\begin{theo}{\cite[Theorem 2.3.10]{clarke1990optimization}}
		\label{The_Clarke}
		For any given $X \in \bb{R}^{n\times p}$, it holds that
		\begin{equation*}
			\partial (f\circ \A)(X) \subseteq \JA(X)[\partial f(\A(X))].
		\end{equation*}
		Moreover, when $\A$ maps any neighborhood of $X$ to a  set which is dense in a neighborhood of  $\A(X)$, we have
		\begin{equation*}
			\partial (f\circ \A)(X) = \JA(X)[\partial f(\A(X))].
		\end{equation*}
	\end{theo}

	Finally, we define the following constants for the theoretical analysis of \ref{NEPen},
	\begin{itemize}
		\item $M_0 := \sup_{X \in \Omega_1  } ~  f(X) - \inf_{X \in \Omega_1 } f(X)$;
		\item $M_1 := \sup_{X \in \Omega_1, W \in \partial f(X) } ~ \norm{W}\ff.$
	\end{itemize}
	It is worth mentioning that both $M_0$ and $M_1$ are independent of the penalty parameter $\beta$.

	\section{Theoretical Properties}

	\subsection{Preliminary lemmas}
	In this subsection, we first present the expression for the Jacobian of the nonlinear mapping $\A$, which could be directly derived from the expression of $\A$. Here we denote $\JA(X)$ as the Jacobian of the mapping $\A$ at $X$, which can be regarded as a linear mapping from $\bb{R}^{n\times p}$ to $\bb{R}^{n\times p}$.

	\begin{prop}
		\label{Prop_ATA}
		For any given $X \in \bb{R}^{n\times p}$, it holds that
		\begin{equation} \label{eq-ATA}
			\A(X)\tp\A(X) - I_p = \frac{1}{64}
			(X\tp X - I_p)^3 \big( 9 (X\tp X)^2 - 33 (X\tp X) + 64 I_p\big).
		\end{equation}
		Moreover, we have for any $X\in \Omega_1$,
		\begin{equation*}
			\norm{\A(X)\tp\A(X) - I_p}_{\mathrm{F}} \,\leq \,
			\norm{X\tp X - I_p}_{\mathrm{F}}^3.
		\end{equation*}
	\end{prop}
	\begin{proof} By using the  economical SVD of $X = U\Sigma V\tp$, we have that
	$$
	\A(X)\tp \A(X) -I_p  = \frac{1}{64} V \Big(
	\Sigma^2 (15 I_p - 10 \Sigma^2 + 3 \Sigma^4)^2 -64 I_p \Big) V^T.
	$$
	By considering the polynomial function $p(x) = x (15-10 x+ 3x^2)^2 - 64$ and noting
	that $p(1) = p'(1) = p''(1) = 0$, we can factorize it as
	$p(x) = (x-1)^3(9x^2 -33 x + 64).$ From here, we can readily obtain the required
	result in \eqref{eq-ATA}.

	Next, for any $X\in \Omega_1$ with the economical SVD of $X = U\Sigma V\tp$, we have
	that $\Sigma^2 \preceq 2 I_p$, and
	\begin{eqnarray*}
	\norm{\A(X)\tp\A(X) - I_p}_{\mathrm{F}} &\leq& \frac{1}{64}\norm{X\tp X - I_p}_{\mathrm{F}}^3
	\norm{ 9 (X\tp X)^2 - 33 (X\tp X) + 64 I_p}_2
	\\
	& =& \frac{1}{64}\norm{X\tp X - I_p}_{\mathrm{F}}^3\; \max_{0\leq x \leq 2}\{
	 |9 x^2 - 33 x + 64|\}
	 \\
	 & = & \norm{X\tp X - I_p}_{\mathrm{F}}^3.
	\end{eqnarray*}
	This completes the proof.
	\end{proof}

	\begin{prop}
		\label{Prop_diff_A}
		For any $X\in \bb{R}^{n\times p}$, it holds that for any $D\in \bb{R}^{n\times p}$
		\begin{equation*}
			\JA(X)[D] = \frac{1}{8}D \left( 15 I_p - 10 X\tp X + 3(X\tp X)^2\right) - X \Phi(X\tp D) + \frac{3}{2}X\Phi(\Phi(X\tp D)(X\tp X - I_p)).
		\end{equation*}
		Moreover, for any $X \in \ca{S}_{n,p}$, we have
		\begin{equation*}
			\JA(X)[D] = D - X\Phi(D\tp X).
		\end{equation*}
	\end{prop}
	The statements in Proposition \ref{Prop_diff_A} can be verified by straightforward calculations, and hence its proof is omitted.

	\begin{lem}
		For any given $X\in \bb{R}^{n\times p}$, the mapping $\JA(X)$ is self-adjoint.
	\end{lem}

	\begin{proof}
		For any $D,W\in \bb{R}^{n\times p}$, it is easy to verify that
		\begin{align*}
			& \tr\left(W\tp D \left( {15} I_p - {10}X\tp X + {3} (X\tp X)^2 \right)\right) = \tr\left(D\tp W \left( {15} I_p - {10}X\tp X + {3} (X\tp X)^2 \right)\right),\\
			& \tr\left( W\tp X \Phi(X\tp D) \right) = \tr\left( D\tp X \Phi(X\tp W) \right).
		\end{align*}
		Moreover, through direct calculation, we achieve
		\begin{equation*}
			\begin{aligned}
				&\tr\left(W\tp X\Phi(\Phi(X\tp D)X\tp X)\right)= \tr\left(\Phi(X\tp W)\Phi(X\tp D)X\tp X\right)\\
				={}& \tr\left(X\tp X\Phi(X\tp W)\Phi(X\tp D)\right) =  \tr\left( D\tp X\Phi(\Phi(X\tp W)X\tp X) \right).
			\end{aligned}
		\end{equation*}
		Therefore, we conclude that
		\begin{equation*}
			\begin{aligned}
				&\inner{\JA(D),W} \\
				={}& \inner{\frac{1}{8}D \left( 15 I_p - 10 X\tp X + 3
				(X\tp X)^2 \right) - \frac{5}{2}X \Phi(X\tp D) +\frac{3}{2}X\Phi(\Phi(X\tp D)X\tp X),W}\\
				={}& \inner{\frac{1}{8}W \left( 15 I_p - 10X\tp X + 3 (X\tp X)^2 \right) - \frac{5}{2}X \Phi(X\tp W) +\frac{3}{2}X\Phi(\Phi(X\tp W)X\tp X),D}\\
				={}& \inner{\JA(W),D},
			\end{aligned}
		\end{equation*}
		and this completes the proof.
	\end{proof}

	\begin{lem}
		\label{Le_dense_surjective}
		$\A$ is a homeomorphism from $\bb{R}^{n\times p}$ to $\bb{R}^{n\times p}$. Moreover, for any given $X \in \bb{R}^{n\times p}$ and any $\varepsilon > 0$, there exists a $\delta > 0$ such that $\ca{B}(\A(X), \delta) \subset \A(\ca{B}(X, \varepsilon) )$.
	\end{lem}
	\begin{proof}
		Consider the auxiliary function $\psi(t) = t(\frac{15}{8} - \frac{5}{4}t^2 + \frac{3}{8}t^4 )$. Since $\psi'(t) = \frac{15}{8}(t^2 - 1)^2$, we can conclude that $\psi(t)$ is non-decreasing in $\bb{R}$ and $\psi'(t)$ is positive except for $t = \pm 1$. As a result, $\psi(t)$ is an injection from $\bb{R}$ to $\bb{R}$, and it is easy to verify that $\psi$ is a bijection and  $\psi^{-1}$ is continuous in $\bb{R}$.

		Let $X = U\Sigma V\tp$ be the rank-revealing singular value decomposition of $X$, where both $U \in \bb{R}^{n\times p}$ and $V \in \bb{R}^{p\times p}$ are orthogonal matrices, and $\Sigma \in \bb{R}^{p\times p}$ is diagonal. Then from the definition of $\psi$, we can  conclude that $\A(X) = U \psi(\Sigma)V\tp$, and $\A$ is an injection. Moreover,  for any given $Y \in \bb{R}^{n\times p}$ and its singular value decomposition $Y = U_Y \Sigma_Y V_Y\tp$, from the definition of $\psi^{-1}$, it holds that $Y = \A(U_Y \psi^{-1}(\Sigma_Y) V_Y\tp)$.  As a result, $\A$ is a bijection and $\A^{-1}$ is well-defined and continuous in $\bb{R}^{n\times p}$.

		Therefore, $\A$ maps any open set to an open set due to the continuity of $\A^{-1}$. As a result, we can conclude that for any $X \in \bb{R}^{n\times p}$ and any $\varepsilon > 0$,  $\A(\ca{B}(X, \varepsilon) )$ is an open set that contains $\A(X)$. Therefore, there exists $\delta > 0$ such that $\ca{B}(\A(X), \delta) \subset \A(\ca{B}(X, \varepsilon) )$ and thus we complete the proof.

	\end{proof}

	\begin{prop}
		\label{Prop_chain_rule}
		For any given $X \in \bb{R}^{n\times p}$,
		\begin{equation*}
			\partial g(X) = \JA(X)[\partial f(\A(X))].
		\end{equation*}
		Moreover, we have
		\begin{equation*}
			\partial h(X) = \JA(X)[\partial f(\A(X))] + \beta X(X\tp X - I_p).
		\end{equation*}
	\end{prop}
	\begin{proof}
		From Lemma \ref{Le_dense_surjective}, and Theorem \ref{The_Clarke} (Theorem 2.3.10 in \citep{clarke1990optimization}), we  have that the chain rule holds for $f(\A(X))$ and
		this leads to the desired result.
	\end{proof}
	Proposition \ref{Prop_chain_rule} illustrates that the Clarke subdifferential of $g$ is easy to achieve from $\partial f$ even if $f$ is only assumed to be locally Lipschitz continuous.

	\begin{rmk}
		From the definition of $h(X)$ and Lemma \ref{Prop_chain_rule}, we can easily verify that for any $X \in \ca{S}_{n,p}$, it holds that
		\begin{equation}
			h(X) = f(X),\quad \text{and }\quad  \partialProj f(X) = \partial h(X).
		\end{equation}
	\end{rmk}

	\subsection{Relationships on stationary points}
	\label{Subsection_equivalence_stationary}
	In this section, we build up the relationships between the first-order stationary points of \ref{Prob_Ori} and \ref{NEPen}.

	\begin{lem}
		For any given $X \in \bb{R}^{n\times p}$ and any $W \in \bb{R}^{n\times p}$, it holds that
		\begin{equation*}
			\inner{J_{\A}(X) [W], X(X\tp X - I_p)} = \frac{15}{8} \inner{\Phi(X\tp W), (X\tp X - I_p)^3}.
		\end{equation*}
	\end{lem}
	\begin{proof}
		From Proposition \ref{Prop_diff_A}, we can obtain after some tedious calculations that
		\begin{equation*}
			\begin{aligned}
				&\inner{J_{\A}(X) [W], X(X\tp X - I_p)}
				=\inner{X\tp J_{\A}(X) [W],(X\tp X - I_p)}\\
				=&\inner{\Phi(X\tp W),\frac{15}{8}((X\tp X)^2-2X\tp X+I_p)(X\tp X - I_p)} \\
				=& \frac{15}{8}\inner{\Phi(X\tp W), (X\tp X - I_p)^3},
			\end{aligned}
		\end{equation*}
		and  the proof is completed.
	\end{proof}

	\begin{theo}
		\label{The_Equivalence_local}

		For any given $\beta > 0$, $r \in \left(0,  \frac{\beta}{2\beta + 8M_1} \right)$, and any $X \in \Omega_{r}$, it holds that
		\begin{equation*}
			\mathrm{dist}(0, \partial h(X)) \geq \frac{\beta}{4} \norm{X\tp X - I_p}\ff.
		\end{equation*}
		Furthermore, when $X^* \in \Omega_{r}$ is a first-order stationary point of \ref{NEPen}, then $X^*$ is a first-order stationary point of \ref{Prob_Ori}.
	\end{theo}
	\begin{proof}
		Suppose $X$ is not feasible, then $\norm{ {X}\tp X - I_p }\ff > 0$.
		Since $X\in \Omega_r$,  it holds that  $\norm{ {X}\tp X - I_p }\ff \leq r < 1$ and
		\begin{equation*}
			\sigma_{\min}({X}\tp X) \geq 1 - \sigma_{\max}({X}\tp X - I_p) \geq 1 - \norm{{X}\tp X - I_p}\ff.
		\end{equation*}
		Then for any $D \in \partial h(X)$, by Proposition \ref{Prop_chain_rule}
		 there exists $W \in \partial f(\A(X))$ such that $D = J_{\A}(X)[W] + \beta  X({X}\tp X - I_p)$. Moreover, by Proposition \ref{Prop_ATA}, $\A(X)\in \Omega_{r^3} \subset \Omega_1$, and
		 hence $\norm{W}\ff \leq M_1$. In addition, since $\norm{X}_2 \leq 2$, we get
		 $\norm{\Phi(X\tp W)}\ff \leq \norm{X}_2 \norm{W}\ff \leq 2 M_1.$
		Now we can conclude that
		\begin{equation*}
			\begin{aligned}
				&2\norm{D}\ff \norm{X\tp X - I_p}\ff \geq \norm{D}\ff \norm{X}_2\norm{X\tp X - I_p}\ff
				\geq \inner{D, {X}({X}\tp {X} - I_p)}\\
				\geq {}&   \beta  \inner{{X}({X}\tp {X} - I_p),{X}({X}\tp {X} - I_p) } - \left| \inner{ J_{\A}(X)[W],  {X}({X}\tp {X} - I_p)} \right| \\
				\geq{}& \beta \sigma_{\min}({X}\tp X ) \norm{{X}\tp X - I_p}\fs
				- \frac{15}{8} \left|\inner{\Phi({X}\tp W), (X\tp X - I_p)^3}\right| \\
				\geq{}& \beta \sigma_{\min}({X}\tp X ) \norm{{X}\tp X - I_p}\fs - \frac{15}{8} \norm{\Phi({X}\tp W)}_2 \norm{X\tp X - I_p}\ff^3\\
				\geq{}& \beta\norm{{X}\tp X - I_p}\fs - \left(\frac{15 M_1}{4} + \beta\right)\norm{{X}\tp X - I_p}\ff^{3}\\
				\geq {}&\beta\norm{{X}\tp X - I_p}\fs \left( 1 - \frac{\beta + 4M_1}{\beta} \norm{{X}\tp X - I_p}\ff \right) \\
				\geq{}& \frac{\beta}{2} \norm{{X}\tp X - I_p}\fs.
			\end{aligned}
		\end{equation*}
		Therefore, it holds that
		\begin{equation*}
			\norm{D}\ff \geq \frac{\beta}{4} \norm{{X}\tp X - I_p}\ff.
		\end{equation*}
		By the arbitrariness of $D \in \partial h(X)$, we get
		\begin{equation*}
			\mathrm{dist}(0, \partial h(X)) \geq \frac{\beta}{4} \norm{X\tp X - I_p}\ff.
		\end{equation*}
		Furthermore, when $X^* \in \Omega_{r}$ is a first-order stationary point of \ref{NEPen}, then
		\begin{equation*}
			0 = \mathrm{dist}(0, \partial h(X^*)) \geq \frac{\beta}{4} \norm{{X^*}\tp {X^*} - I_p}\ff \geq 0,
		\end{equation*}
		which implies that  ${X^*}\tp X^* = I_p$ and this completes the proof.
	\end{proof}

	\subsection{Estimating stationarity}
	Directly applying unconstrained optimization approaches to solve \ref{NEPen} usually yields an infeasible sequence and terminates at an infeasible point. However, sometimes we expect mild accuracy in the substationarity, but pursue high accuracy in the feasibility in minimizing \ref{NEPen}. To this end, we impose an orthonormalization post-process after obtaining the solution $X$  by applying unconstrained approaches to solve \ref{NEPen}. Namely,
	\begin{eqnarray}\label{orth}
		\Xo := \ProjS(X),
	\end{eqnarray}
	where $\ProjS: \bb{R}^{n\times p} \to \ca{S}_{n,p}$ is the projection on Stiefel manifold defined in Section 1.4.

	\begin{lem} \label{lem-proj}
		For any given $X \in \bb{R}^{n\times p}$, it holds that
		\begin{equation*}
			\norm{X - \ProjS(X)}\ff \leq \norm{X\tp X - I_p}\ff.
		\end{equation*}
	\end{lem}
	\begin{proof}
		By the economical SVD of $X$, that is, $X = U \Sigma V\tp$ with
		$\Sigma = {\rm diag}(\sigma)$, we can achieve that
		\begin{equation*}
			\norm{X - \ProjS(X)}\ff = \norm{\Sigma - I_p}\ff \leq \norm{\Sigma^2 - I_p}\ff = \norm{X\tp X - I_p}\ff.
		\end{equation*}
		The inequality above holds simply because $|x-1| \leq |x+1||x-1|$
		for any $x\geq 0$.
	\end{proof}

	\begin{lem} \label{lem-A-proj}
		For any given $X \in \Omega_{1/2}$, it holds that
		\begin{equation} \label{eq-A-proj}
			\norm{\A(X) - \ProjS(X)}\ff \leq 4\norm{X\tp X - I_p}\ff^3.
		\end{equation}

		Moreover
		\begin{eqnarray} \label{eq-fA-fproj}
		 | f(\A(X)) - f(\ProjS(X))| \;\leq \; M_1 \norm{\A(X) - \ProjS(X)}\ff.
		\end{eqnarray}
	\end{lem}
	\begin{proof}
		By the economical SVD of $X$, that is, $X = U \Sigma V\tp$, we obtain that
		$$
		  \A(X) - \ProjS(X) = \frac{1}{8} U\Big( 3 \Sigma^5 - 10 \Sigma^3 + 15 \Sigma - 8 I_p
		  \Big) V^T.
		$$
		Note that $X\in \Omega_{1/2}$ implies that $0\prec \Sigma \preceq \frac{3}{2} I_p$.
		Next, we have that
		\begin{equation*}
			\begin{aligned}
				&\norm{\A(X) - \ProjS(X)}\ff = \frac{1}{8}\norm{3\Sigma^5 - 10 \Sigma^3
				+ 15\Sigma - 8 I_p }\ff \\
				={}& \frac{1}{8}\norm{\left( \Sigma - I_p \right)^3 \left( 3 \Sigma^2 + 9\Sigma + 8I_p   \right)}\ff \\
				\leq {}& \frac{1}{8}
				\norm{ 3 \Sigma^2 + 9\Sigma + 8I_p   }_2 \norm{\left( \Sigma - I_p \right)^3}\ff
				 \\
				\leq {}&  \left(\frac{1}{8} \max_{0\leq x \leq 3/2} \{
				 |3 x^2 + 9 x + 8|\}\right) \norm{\Sigma - I_p}\ff^3
				\leq 4\norm{\Sigma^2 - I_p}\ff^3 = 4\norm{X\tp X - I_p}\ff^3.
			\end{aligned}
		\end{equation*}
		Note that we used Lemma \ref{lem-proj} in the last inequality.

		Next we prove \eqref{eq-fA-fproj}.
		Since $f$ is locally Lipschitz continuous, by \cite[Theorem 2.3.7]{clarke1990optimization},
		there exists $t\in (0,1)$ such that
		$f(\A(X)) - f(\ProjS(X)) \in \{ \langle W, \, \A(X) - \ProjS(X)\rangle : W \in \partial f(Z_t)\}$, where $Z_t =   t \A(X) + (1-t) \ProjS(X)$. Thus
		$$
		 |f(\A(X)) - f(\ProjS(X))| \leq  \norm{\A(X) - \ProjS(X)}\ff \sup_{W\in \partial f(Z_t)} \norm{W}\ff
		$$
		By some tedious calculations and using the result in \eqref{eq-A-proj}, we can show that
		$Z_t \in \Omega_1$, and hence $\sup_{W\in \partial f(Z_t)} \norm{W}\ff \leq M_1$.
		This completes the proof.

	\end{proof}

	The following proposition illustrates that the post-processing procedure \eqref{orth}  can further reduce the function value simultaneously if the current point is close to the Stiefel manifold.

	\begin{prop}
		\label{Le_ortho_fval}

		For any given
		$X \in \Omega_{1/2}$,  it holds that
		\begin{equation*}
			h(\ProjS(X)) \leq h(X) - \left( \frac{\beta}{4} -  4M_1\norm{X\tp X - I_p}\ff \right)\norm{X\tp X - I_p}\fs.
		\end{equation*}
	\end{prop}
	\begin{proof}
		\begin{equation*}
			\begin{aligned}
				&h(X) - h(\ProjS(X)) = f(\A(X)) - f(\ProjS(X)) +  \frac{\beta}{4}\norm{X\tp X - I_p}\fs \\
				\overset{(i)}{\geq}{}& -M_1 \norm{\A(X) - \ProjS(X)}\ff + \frac{\beta}{4}\norm{X\tp X - I_p}\fs  \\
				\geq{}& \left( \frac{\beta}{4} -  4M_1 \norm{X\tp X - I_p}\ff \right)\norm{X\tp X - I_p}\fs.
			\end{aligned}
		\end{equation*}
		Here $(i)$ follows from \eqref{eq-fA-fproj} in Lemma \ref{lem-A-proj}.
	\end{proof}

	\begin{theo}
		\label{The_Equivalence_local_minimizer}

		For any given $X^* \in \Omega_{1/2}$, suppose  $X^*$ is a local minimizer of \ref{NEPen}, then  $X^*$ is a local minimizer of \ref{Prob_Ori}. Moreover, when $\beta \geq 16M_1$, any local minimizer $X^*$ of \ref{Prob_Ori} is a local minimizer of \ref{NEPen}.
	\end{theo}
	\begin{proof}
		Suppose $X^* \in \Omega_{1/2}$ is a local minimizer of \ref{NEPen}. Then $0 \in \partial h(X^*)$ and Theorem \ref{The_Equivalence_local} further implies that $X^* \in \ca{S}_{n,p}$. Notice that $h(X) = f(X)$ holds for any $X \in \ca{S}_{n,p}$, we immediately obtain that $X^*$ is a local minimizer of \ref{Prob_Ori}.

		On the other hand, when $X^* \in \ca{S}$ is a local minimizer of \ref{Prob_Ori},  then there exists a $\gamma \in (0, 1)$ such that $f(Z) \geq f(X^*)$ holds for any $Z \in \ca{S}_{n,p}$ satisfying $\norm{Z - X^*}\ff \leq \gamma$. Then for any $Y \in \bb{R}^{n\times p}$ such that  $\norm{Y - X^*}\ff \leq \frac{\gamma}{2}$ and $\norm{Y\tp Y - I_p}\ff \leq \frac{\gamma}{2}$, we first have that
		\begin{equation*}
			\norm{\ProjS(Y) - X^*}\ff \leq \norm{\ProjS(Y) - Y}\ff + \norm{Y-X^*}\ff \leq \frac{\gamma}{2} + \frac{\gamma}{2} =  \gamma.
		\end{equation*}
		Then from Proposition \ref{Le_ortho_fval} we can conclude that
		\begin{equation*}
			\begin{aligned}
				h(Y) - h(X^*) ={}& h(Y) - h(\ProjS(Y)) + h(\ProjS(Y)) - h(X^*) \\
				\geq{}& \left( \frac{\beta}{4} -  4M_1 \norm{Y\tp Y - I_p}\ff \right) \norm{Y\tp Y - I_p}\fs \geq 0,
			\end{aligned}
		\end{equation*}
		and complete the proof.
	\end{proof}

	\begin{rmk}
		For the following optimization problem that minimizes a nonsmooth objective function over product of multiple Stiefel manifolds and Euclidean spaces,
			\begin{equation}
				\label{Eq_Prob_Multi_Ori}
				\begin{aligned}
					\min_{Y_i \in \bb{R}^{n_i\times p_i}, Z \in \bb{R}^l} \quad& f(Y_1,...,Y_N, Z)\\
					\text{s. t.} \quad & Y_i\tp Y_i = I_{p_i}, ~\text{for~} i = 1,...,N,
				\end{aligned}
			\end{equation}
			the corresponding constraint dissolving function can be formulated as
			\begin{equation}
				\label{Eq_Prob_Multi_Pen}
				f(\A_1(Y_1),...,\A_N(Y_N), Z) + \frac{\beta}{4} \left( \sum_{i = 1}^N \norm{Y_i\tp Y_i - I_{p_i}}\fs \right),
			\end{equation}
			where $\A_i(Y_i) := \frac{1}{8} Y_i \left( 15 I_{p_i} - 10 Y_i\tp  Y_i + 3(Y_i\tp Y_i)^2 \right)$.
			The relationship between \eqref{Eq_Prob_Multi_Ori} and \eqref{Eq_Prob_Multi_Pen}  can be established in the same way as in the proofs in Section 3.

		Moreover, the exactness of \ref{NEPen} can be further extended to optimization problems over the generalized Stiefel manifold,
		\begin{equation}
			\label{Prob_Gen}
			\begin{aligned}
				\min_{X \in \bb{R}^{n\times p}} \quad& f(X)\\
				\text{s. t.} \quad & X\tp B X = I_p,
			\end{aligned}
		\end{equation}
		where $B \in \bb{R}^{n\times n}$ is a symmetric positive definite matrix.  For \eqref{Prob_Gen}, we can consider the following exact penalty function $\hsharp$,
		\begin{equation}
			\label{penalty_function_general}
			\hsharp(X) := f(\A^\sharp(X)) + \frac{\beta}{4} \norm{X\tp BX - I_p}\fs,
		\end{equation}
		where
		\begin{equation}
			\A^\sharp(X):= \frac{1}{8}X \left( 15 I_p - 10 X\tp B X + 3 (X\tp BX)^2 \right).
		\end{equation}
		Similar to the proof of Lemma \ref{Le_dense_surjective}, it is easy to verify that $\A^{\sharp}$ is a homeomorphism from $\bb{R}^{n\times p}$ to $\bb{R}^{n\times p}$. Therefore, we can the prove the exactness of \eqref{penalty_function_general} together with the results in Section 2 and Section 3 through the same approach as the proofs in this paper.
	\end{rmk}

	\section{Subgradient Methods}

		In this section, we discuss how to develop subgradient methods for \ref{Prob_Ori} and establish its convergence properties. Subgradient method and its variants play important roles in minimizing nonsmooth functions that are not regular, in particular, not prox-friendly.  Recently, \cite{davis2020stochastic} shows the global convergence for subgradient methods in minimizing those functions that {\it admit a chain rule}, which is defined later in Definition \ref{Defin_Chain_rule}. Such properties are essential for the so-called ``descent condition'' \cite[Assumption D]{davis2020stochastic} for subgradient methods, hence are fundamental for establishing their convergence properties for solving nonconvex nonsmooth optimization, as illustrated in various existing works \cite{davis2020stochastic,bolte2021conservative,castera2021inertial}. Without assuming such property, as far as we know, there is no technique to analyze these theoretical properties of subgradient methods in such general situations. Therefore, we first introduce the concept of functions that admit a chain rule in Section \ref{Subsection_Whitney}. Moreover,  we show that functions satisfying such property can cover all the problems considered in this paper, and remarkably \ref{NEPen} admits a chain rule whenever the objective function $f$ admits a chain rule. Furthermore,  we prove several important preliminary lemmas in Section \ref{Subsection_4_Lemmas}, and  propose a unified framework of subgradient methods in Section \ref{Subsection_4_subgradient_method}.

	\begin{rmk}
		By transforming \ref{Prob_Ori} into \ref{NEPen}, various existing approaches in unconstrained nonsmooth nonconvex optimization can be directly implemented to solve \ref{Prob_Ori}. These approaches include gradient sampling methods  \cite{burke2005robust,kiwiel2007Convergence,davis2021gradient}, normalized subgradient methods \cite{zhang2020complexity}, nonsmooth second-order approaches \cite{milzarek2019stochastic,asl2021analysis,yang2021stochastic}, etc. Furthermore, the convergence properties of those approaches, including the global convergence and iteration complexity, directly follow the  existing related works \cite{davis2019stochastic,zhang2020complexity,davis2021proximal,davis2021gradient,gratton2020algorithm}.
	\end{rmk}

	\subsection{Preliminaries}
	\label{Subsection_Whitney}

	In this subsection, we first present the concept of the functions that admit a chain rule. Then we introduce the concept of Whitney stratifiable functions, which plays an important role in nonsmooth analysis and optimization \cite{bolte2007clarke}. We show that all the Whitney stratifiable functions and Clarke regular functions admit a chain rule. Furthermore, we provide the definition for several important classes of functions, including semi-algebraic functions, semi-analytic functions, and functions that are definable in an $o$-minimal structure, all of which are contained in the class of Whitney stratifiable functions.

	\begin{defin}
		An absolutely continuous curve is a continuous $\gamma: [0,1] \to \bb{R}^{n\times p}$ whose derivative exists almost everywhere in $\bb{R}$, and $\gamma(t) - \gamma(0)$ is Lebesgue integral of $\gamma'$ between $0$ and $t$ for all $t \in \bb{R}$, i.e.,
		\begin{equation}
			\gamma(t) = \gamma(0) + \int_{0}^t \gamma'(\tau) \mathrm{d} \tau \qquad \text{for all $t \in \bb{R}$}.
		\end{equation}
	\end{defin}

	\begin{defin}[Deinfition 5.1 in \cite{davis2020stochastic}]
		\label{Defin_Chain_rule}
		We say a locally Lipschitz continuous function $f$ admits a chain rule if for any absolutely continuous curve $\gamma: \bb{R}_+ \to \bb{R}^{n\times p}$, the equality
		\begin{equation}
			(f\circ \gamma)'(t) = \inner{\partial f(\gamma(t)), \gamma'(t)}
		\end{equation}
		holds for a.e. $t \geq 0$.
	\end{defin}

	\begin{lem}[\cite{bolte2021conservative}]
		Suppose $f$ admits a chain rule, then for any absolutely continuous curve $\gamma: [0,1] \to \bb{R}^{n\times p}$, it holds that
		\begin{equation}
			f(\gamma(1)) - f(\gamma(0)) = \int_{0}^{1} \max_{V \in \partial f(\gamma(t))} \inner{V, \gamma'(t)} \mathrm{d} t = \int_{0}^{1} \min_{V \in \partial f(\gamma(t))} \inner{V, \gamma'(t)} \mathrm{d} t.
		\end{equation}
	\end{lem}

	The authors in \cite{davis2020stochastic,bolte2021conservative} prove the global convergence for applying a subgradient method to minimize locally Lipschitz continuous function that admits a chain rule. And it is also proved  that Clarke regular functions admit a chain rule\cite{davis2020stochastic}.
		\begin{prop}[Lemma 5.4 in \cite{davis2020stochastic}]
			Any locally Lipschitz function that is Clarke regular admits a chain rule.
		\end{prop}

		In the rest of this subsection, we focus on another broad class of functions  known as Whitney stratifiable functions.
	Before giving a formal definition, let us fix some notations. A set $M \subset \bb{R}^{n\times p}$ is a $C^s$-smooth manifold if there is an integer $r \in \bb{N}$ such that around any point
	$X \in M$, there is a neighborhood $U$ and a $C^s$-smooth mapping $F: U \mapsto \bb{R}^{np-r}$ such that the Jacobian of $F$ is full rank at any $ Y \in U$ and $M \cap U= \{ V \in U: F(V) = 0 \}$. And we denote the tangent space and normal space of $M$ at $X$ as $\ca{T}_{M}(X)$ and $\ca{N}_{M}(X)$, respectively.

	\begin{defin}{\cite[Definition 5.6]{davis2020stochastic}}

		A Whitney $\ca{C}^s$-stratification $\ca{E}$ of a set $Q \subset \bb{R}^{n\times p}$ is a partition of $Q$ into finitely many nonempty $C^s$ manifolds, called strata, satisfying the following compatibility conditions.
		\begin{itemize}
			\item {\bf Frontier condition:} For any two strata $L$ and $M$, the following implication holds
			\begin{equation*}
				L \cap \mathrm{cl}(M) \neq \emptyset \quad \Longrightarrow \quad L \subset \mathrm{cl}(M).
			\end{equation*}
			\item {\bf Whitney condition:} For any sequence of points $\{Z_k\}$ in stratum $M$ converging to a point $\bar{Z}$ in stratum $L$, if the corresponding normal vectors $\{V_k\} $ chosen by $ V_k \in \ca{N}_{M}(Z_k)$ converge to a vector $V$, then $V \in \ca{N}_{L}(\bar{Z})$.
		\end{itemize}
	\end{defin}

	Moreover, we say that a function $f: \bb{R}^{n\times p} \to \bb{R}$ is Whitney $\ca{C}^s$-stratifiable if its graph admits a Whitney $\ca{C}^s$-stratification.  The following theorem illustrates that Whitney stratification implies the chain rule for locally Lipschitz continuous functions.
	\begin{theo}[Theorem 5.8 in \cite{davis2020stochastic}]
		Any locally Lipschitz function $f:\bb{R}^{n\times p} \to \bb{R}$ that is Whitney $C^1$-stratifiable admits the chain rule.
	\end{theo}

	It is worth mentioning that  verifying the path differentiability for nonsmooth functions is usually an easy task.
	As shown in \cite{davis2020stochastic}, any Clarke regular function admits a chain rule. Moreover, Whitney $\ca{C}^1$-stratifiable functions contain a broad class of functions, including semi-algebraic and semi-analytic functions. 	A closed set $Q$ is called semi-analytic if it can be written as a finite union of sets, each takes the form as
	\begin{equation}
		\{ x \in \bb{R}^{n\times p}: p_i(x) \leq 0\quad \text{for $i = 1,...l$} \},
	\end{equation}
	for some real-analytic functions $p_1,...,p_l$ defined on $ \bb{R}^{n\times p}$. If all the function $p_i$ are polynomials, then we say that $Q$ is semi-algebraic. As mentioned in \cite{law1965ensembles}, any semi-analytic set admits a Whitney $C^{\infty}$ stratification.

	Another important class of functions is characterized based on the concept of o-minimal structure.   We first present the definition for o-minimal structure, which follows the definitions in \cite{van1996geometric,davis2020stochastic}:
	\begin{defin}
		An o-minimal structure is a collection of set $\ca{O} = \{\ca{O}_k\}_{k = 1,2,...}$ where each $\ca{O}_k$ is a family of subsets of $\bb{R}^d$ such that for each $k \geq 0$:
		\begin{enumerate}
			\item $\ca{O}_k$ is closed under complementation, finite union, finite intersection and contains $\bb{R}^k$.
			\item If $A$ belongs to $\ca{O}_k$, then $A \times \bb{R}$ and $\bb{R} \times A$ belongs to $\ca{O}_{k+1}$.
			\item If $\pi: \bb{R}^k \times \bb{R} \to \bb{R}^k$ denotes the coordinate projection onto $\bb{R}^k$, then for any $A$ in $\ca{O}_{k+1}$, the set $\pi(A)$ belongs to $\ca{O}_k$.
			\item $\ca{O}_k$ contains the family of real algebraic subsets of $\bb{R}^k$, that is, the set of the form $\{ x \in \bb{R}^k: p(x) = 0 \}$ for a polynomial $p(x)$ on $\bb{R}^k$.
			\item The elements of $\ca{O}_1$ are exactly the finite unions of intervals.
		\end{enumerate}

		Moreover, the sets $A$ belonging to $\ca{O}_d$ for some $d \geq 1$ are called definable in the o-minimal structure.
	\end{defin}

		\begin{defin}
			A locally Lipschitz function $f$ is said to be definable if its graph is definable in an $o$-minimal structure.
		\end{defin}

		As illustrated in \cite{van1996geometric}, any function definable in an $o$-minimal structure admits a Whitney $\ca{C}^p$- stratification for any $p \geq 1$, and hence admits a chain rule. It  straightforwardly holds from Tarski–Seidenberg theorem \cite{bierstone1988semianalytic} that any semi-algebraic function is definable.
		Moreover, beyond semialgebraicity, \cite{wilkie1996model} shows that there is an o-minimal structure that simultaneously contains both the
		graph of the exponential function and all semi-algebraic sets. As a result, various common activation functions and loss functions, such as sigmoid, softplus, RELU, $\ell_1$-loss, MSE loss, hinge loss, logistic loss and cross-entropy loss, are all definable. Additionally, as shown in \cite{davis2020stochastic,bolte2021conservative}, any composition of two definable functions is definable, which further illustrates that the loss function of a neural network built from definable activation functions is definable and thus admits a chain rule.

	Therefore, the objective functions for the applications mentioned in \cite{absil2017collection,chen2018proximal,lerman2018overview,tsakiris2018dual,chen2019manifold,xiao2020l21,xiao2021penalty} are  Whitney $\ca{C}^1$-stratifiable. As a result, we can develop subgradient methods to minimize those functions and prove their convergence properties based on the framework presented in \cite{davis2020stochastic}.

	\begin{assumpt}
		\label{Assumption_f}
		For the objective function $f$ in \ref{Prob_Ori}, we assume,
		\begin{enumerate}
			\item The objective function $f$ in \ref{Prob_Ori} admits a chain rule;
			\item The set $\{ f(X) : 0 \in \partialProj f(X),~ X \in \ca{S}_{n,p} \}$ has empty interior in $\bb{R}$.
		\end{enumerate}

	\end{assumpt}
	It is worth mentioning that the mapping $\A$ is semi-algebraic, then it follows from \cite{bolte2021conservative} that \ref{NEPen} admits a chain rule whenever the objective function $f$ in \ref{Prob_Ori} admits a chain rule.

	\subsection{Preliminary lemmas}
	\label{Subsection_4_Lemmas}

We can conclude from Section \ref{Subsection_equivalence_stationary} that when an unconstrained optimization approach is applied to minimize \ref{NEPen}, if the generated sequence $\{\Xk\}$ is contained in $\Omega_{1/6}$ and converges to some stationary points of \ref{NEPen}, then $\{\Xk\}$ converges to some stationary points of \ref{Prob_Ori}.

	In this subsection, we present two preliminary lemmas to illustrate that keeping the generated sequence restricted in $\Omega_{1/6}$ only requires mild conditions. The following lemma shows that when starting from an initial point in $\Omega_{1/12}$, a specific algorithm yields a sequence $\{ \Xk \}$ such that $h(\Xk) \leq h(X_0)$ holds for any $k \geq 0$, then we can conclude that the sequence $\{\Xk\}$ is restricted in $\Omega_{1/6}$.

	\begin{lem}
			\label{Le_level_bounded}
			For any $Y \in \Omega_{1/12}$ and  $Z \in \Omega_{1} \setminus \Omega_{1/6}$, suppose $\beta \geq 192M_0$, then it holds that  $h(Y) \leq h(Z)$.
		\end{lem}
		\begin{proof}
			For any $Y \in \Omega_{1/12}$, $Z \in \Omega_{1} \setminus \Omega_{1/6}$, we have
			\begin{equation}
				h(Y) - h(Z) \leq \sup_{W \in \Omega_{1} } f(W) - \inf_{W \in \Omega_{1} } f(W) + \frac{\beta}{576}- \frac{\beta}{144} \leq M_0 -\frac{\beta}{192} \leq 0,
			\end{equation}
			and complete the proof.
	\end{proof}
	Lemma \ref{Le_level_bounded} provides  theoretical guarantees for a large number of algorithms,  such as the gradient sampling methods \cite{burke2005robust,kiwiel2007Convergence,davis2021gradient}, the nonsmooth second-order approaches \cite{milzarek2019stochastic,yang2021stochastic}, etc.

	Moreover, when applying subgradient-based algorithms to minimize the penalty function, we usually need to choose $\Dk$ to approximate $\partial g(\Xk)$  at every iterate $\Xk$. The direction $D_k$ is directly chosen from $\partial g(\Xk)$ in deterministic subgradient methods, or with noise in stochastic settings. In the following lemmas, we show that under mild conditions, the generated sequence $\{\Xk\}$ is restricted in $\Omega_{1/6}$. Therefore, the convergence properties of deterministic subgradient methods could be established without assuming that $h$ is bounded below in $\bb{R}^{n\times p}$.

	\begin{lem}
		\label{Le_iter_bounded_l2}
		Suppose $Y_0$ satisfies $Y_0 \in \Omega_{1/6}$ and $\Yk$ is generated by
		\begin{equation*}
			\Ykp = \Yk - \etak (\Dk + \beta \Yk (\Yk \tp \Yk - I_p)),
		\end{equation*}
		where $\Dk \in \bb{R}^{n\times p}$, and there exists a constant $M_2$ such that
		$\norm{\Dk}\ff \leq M_2$
		holds for any $k \geq 0$. Then when $\beta \geq 60M_2$ and $\etak \leq \frac{1}{2\beta}$, $\Yk \in \Omega_{1/6}$ holds for any $k \geq 0$.
	\end{lem}
	\begin{proof}
	Note that when $\norm{\Yk\tp \Yk - I_p}\ff \leq \frac{1}{6}$, we have that $5/6I_p \preceq \Yk\tp \Yk \preceq 7/6 I_p$, and hence $\norm{\Yk}_2 \leq 7/6.$
		By the definition of $\Yk$, let $\hat{D}_k := \Dk + \beta \Yk (\Yk \tp \Yk - I_p)$.
	     We have that
		\begin{equation*}
			\begin{aligned}
				&\norm{\Ykp\tp \Ykp - I_p}\ff
				= \norm{\Yk\tp \Yk - I_p - \etak \hat{D}_k\tp \Yk - \etak  \Yk\tp \hat{D}_k
				+ \etak^2 \hat{D}_k\tp\hat{D}_k  }\ff \\
				\leq{}&
				\norm{\Yk\tp \Yk - I_p - 2\etak \beta \Yk\tp \Yk(\Yk\tp \Yk - I_p) + \etak^2 \beta^2 \Yk\tp\Yk (\Yk \tp \Yk - I_p)^2}\ff \\
				&+ \etak\norm{ -\Dk\tp \Yk - \Yk\tp \Dk + \etak\beta  (\Yk\tp \Yk - I_p) \Yk\tp\Dk
				 + \etak\beta \Dk \tp \Yk(\Yk\tp \Yk - I_p) + \etak \Dk\tp\Dk }\ff\\
				\leq{}&  \norm{ (\Yk\tp \Yk - I_p)\Big((I_p - \etak\beta\Yk\tp \Yk)^2 -\etak^2
				\beta^2\Yk\tp\Yk\Big) }\ff  \\
			   &+ \etak\Big(2\norm{\Dk}\ff\norm{\Yk}_2 + 2\etak\beta\norm{\Yk\tp \Yk - I_p}\ff
				\norm{\Dk}\ff\norm{\Yk}_2 + \etak \norm{\Dk}^2\ff
				\Big)  \\
		\leq{}&		\norm{ (\Yk\tp \Yk - I_p)}\ff \norm{ (I_p - \etak\beta\Yk\tp \Yk)^2 -\etak^2
				\beta^2\Yk\tp\Yk }_2
				+ \etak\Big( 3M_2  + \etak M_2^2\Big) \\
				\leq{}& \left(1 - \frac{5}{6}\etak\beta \right) \norm{\Yk\tp \Yk - I_p}\ff + 3\etak M_2 + \etak^2 M_2^2.
			\end{aligned}
		\end{equation*}
		As a result, when $\norm{\Yk\tp \Yk -I_p}\ff \geq \frac{1}{12}$,
		\begin{equation*}
			\norm{\Ykp\tp \Ykp - I_p}\ff - \norm{\Yk\tp \Yk - I_p}\ff \leq -\frac{5\etak\beta}{72} +  3\etak M_2 + \etak^2 M_2^2 \leq 0.
		\end{equation*}
		Moreover, when $\norm{\Yk\tp \Yk -I_p}\ff \leq \frac{1}{12}$,
		\begin{equation*}
			\norm{\Ykp\tp \Ykp - I_p}\ff \leq \norm{\Yk\tp \Yk - I_p}\ff + 3\etak M_2 + \etak^2 M_2^2 \leq \frac{1}{6}.
		\end{equation*}
		Then, we have  $\norm{\Ykp\tp \Ykp - I_p}\ff \leq \frac{1}{6}$. Therefore, $\norm{\Yk\tp \Yk - I_p}\ff \leq \frac{1}{6}$ holds for $k = 0,1,2,...$ by induction.
	\end{proof}
	Lemma \ref{Le_iter_bounded_l2} illustrates that when we choose the initial point in $\Omega_{1/6}$, the sequence generated by a class of subgradient-based methods are uniformly restricted in $\Omega_{1/6}$.

	\subsection{Convergence properties of subgradient methods}
	\label{Subsection_4_subgradient_method}
	In this subsection, we consider developing subgradient methods for solving \ref{Prob_Ori} by applying existing Euclidean subgradient methods to solve \ref{NEPen}.
	We first present a unified framework (Algorithm \ref{Alg:subgradient_sto}) based on the framework proposed by \cite{davis2020stochastic}.
	Then  we present a basic framework for subgradient methods and its stochastic version.
	Based on the framework presented in Algorithm \ref{Alg:subgradient_sto}, we analyze their convergence properties and show that these subgradient methods are globally convergent.

	\begin{algorithm}[htbp]
		\begin{algorithmic}[1]
			\Require functions $f$ and penalty parameter $\beta$;
			\State Choose an initial guess $X_0 \in \Omega_{1/6}$, set $k=0$;
			\While{not terminate}
			\State Compute $D_k $ as an approximated evaluation for $ \JA(\Xk)[\partial f(\Xk)]$ ;
			\State $\Xkp = \Xk - \eta_k\left( D_k + \beta\Xk (\Xk\tp \Xk - I_p)\right) $;
			\State $k = k+1$;
			\EndWhile
			\State Return $X_k$.
		\end{algorithmic}
		\caption{A framework of subgradient methods for solving \ref{Prob_Ori}.}
		\label{Alg:subgradient_sto}
	\end{algorithm}

	\begin{cond}

		\label{Assumption_Subgradient_sto}
		In Algorithm \ref{Alg:subgradient_sto}, we assume
		\begin{enumerate}

			\item $D_k = G_k + E_k$, where $G_k \in \JA(\Xk)[\partial f(\Xk)]$ and $\lim\limits_{N \to +\infty} \sum\limits_{k = 0}^{N} E_k \eta_k$ exists.
			\item There exists a constant $\tilde{M}$ such that $\sup_{k \geq 0} \norm{D_k}\ff \leq \tilde{M}$, and we choose $\beta \geq \max\{ 16M_1, 60\tilde{M} \}$.
			\item  The sequence $\{\eta_k\}$ satisfies
			\begin{equation*}
				\eta_k > 0, \quad \sum_{k = 0}^{+\infty} \eta_k = \infty, \quad \sum_{k = 0}^{+\infty} \eta_k^2 < +\infty, \quad \sup_{k \geq 0} ~\etak \leq \frac{1}{2\beta}.
			\end{equation*}
		\end{enumerate}
	\end{cond}

	\begin{theo}
		\label{Theo_subgradient_convergence_base}
		Suppose that Assumption \ref{Assumption_f} holds, and Algorithm \ref{Alg:subgradient_sto} satisfies Condition \ref{Assumption_Subgradient_sto}. Then $\{h(\Xk)\}$ converges and every limit point $X^*$ of $\{\Xk\}$ is a stationary point of \ref{Prob_Ori}.
	\end{theo}
	\begin{proof}
		Firstly, Lemma \ref{Le_iter_bounded_l2} illustrates that when $\sup_{k \geq 0} \norm{D_k}\ff \leq \tilde{M}$, $\beta \geq \max\{ 16M_1, 60\tilde{M} \}$, and $\etak \leq \frac{1}{2\beta}$ for any $k \geq 0$ in Algorithm \ref{Alg:subgradient_sto}, it holds that $\{\Xk\}$ is restricted in $\Omega_{1/6}$. Then every cluster point of $\{\Xk\}$ lies in $\Omega_{1/6}$.

		Now we check the validity of Assumption A and Assumption B in \cite{davis2020stochastic}. The fact that $\{\Xk\}$ is restricted in $\Omega_{1/6}$ directly shows that Assumption A(1) and A(2) hold. The Assumption A(3) and A(4) directly follow Condition \ref{Assumption_Subgradient_sto}(2) and \ref{Assumption_Subgradient_sto}(3), respectively. In addition, Assumption A.5 follows quickly from the fact that $\partial h$ is outer-semicontinuous as mentioned in \cite[Lemma 4.1]{davis2020stochastic}.

		Furthermore, Assumption B(1) in \cite{davis2020stochastic} is guaranteed by Assumption \ref{Assumption_f}(2) and Theorem \ref{The_Equivalence_local}, while Assumption B(2) directly follows Assumption \ref{Assumption_f}(1) and \cite[Lemma 5.2]{davis2020stochastic}.

		Therefore, from \cite[Theorem 3.2]{davis2020stochastic} it holds that any cluster point of  $\{\Xk\}$ lies in the set $\{ X \in \Omega: 0 \in \partial h(X) \}$ and $h(\Xk)$ converges. Then it follows from Theorem \ref{The_Equivalence_local} that any point in $\{ X \in \Omega: 0 \in \partial h(X) \}$ is a stationary point of \ref{Prob_Ori}.  Hence we complete the proof.
	\end{proof}

	It is worth mentioning that we can apply Algorithm \ref{Alg:subgradient_sto} to analyze the convergence for stochastic subgradient methods for solving \ref{Prob_Ori}, where the stochasticity is modeled by the sequence $\{E_k\}$ in Condition \ref{Assumption_Subgradient_sto}(2). We follow the standard assumptions that conditioned on the past, each random
	variable $\{E_k\}$ has mean zero and bounded covariance. We present detailed descriptions in Condition \ref{Assumption_Subgradient_sto_new}.

	\begin{cond}
		\label{Assumption_Subgradient_sto_new}
		In Algorithm \ref{Alg:subgradient_sto}, we assume
		\begin{enumerate}
			\item $D_k = G_k + E_k$, where $G_k \in \JA(\Xk)[\partial f(\Xk)]$ and $\{E_k\}$ is a martingale difference sequence with respect to the increasing $\sigma$-fields  $\ca{F}_k := \sigma(X_k, G_j, E_j: j<k)$. That is,  there exists a constant $\hat{M}$, and for  each $k \geq 0$, $E_k$ is measurable with respect to $\ca{F}_k$ and satisfies
				\begin{equation}
					\bb{E}[E_k| \ca{F}_k] = 0, \quad \bb{E}[\norm{E_k}\fs| \ca{F}_k] < \hat{M}.
			\end{equation}
			\item There exists a constant $\tilde{M}$ such that $\sup_{k \geq 0} \norm{D_k}\ff \leq \tilde{M}$ holds almost surely. Moreover, we choose $\beta \geq \max\{ 16M_1, 60\tilde{M} \}$.
			\item  The deterministic sequence $\{\eta_k\}$ satisfies
			\begin{equation}
				\eta_k > 0, \quad \sum_{k = 0}^{+\infty} \eta_k = \infty, \quad \sum_{k = 0}^{+\infty} \eta_k^2 < +\infty, \quad \sup_{k \geq 0} ~\etak \leq \frac{1}{2\beta}.
			\end{equation}
		\end{enumerate}
	\end{cond}

	Condition \ref{Assumption_Subgradient_sto_new}(2) guarantees that sequence $\sum_{k = 0}^{N} \etak E_k$ converges almost surely as mentioned in \cite{davis2020stochastic}. Therefore, Condition \ref{Assumption_Subgradient_sto_new} directly implies Condition \ref{Assumption_Subgradient_sto}, which yields the following corollary on the convergence of the proposed stochastic gradient method.
	\begin{coro}
		Suppose that Assumption \ref{Assumption_f} holds, and Algorithm \ref{Alg:subgradient_sto} satisfies Condition \ref{Assumption_Subgradient_sto_new}. Then with probability $1$,  $\{h(\Xk)\}$ converges  and every limit point $X^*$ of $\{\Xk\}$ is a stationary point of \ref{Prob_Ori}.
	\end{coro}

	\section{Numerical Experiments}

	In this section, we apply \ref{Prob_Ori} to train an orthogonally constrained convolutional neural network (OCNN) and compare its numerical performance with existing practical approaches.  All the numerical experiments in this section are run on a server with Intel Xeon Silver 4110 CPU @ 2.10GHz CPU and NVIDIA Tesla P40 GPU running Pytorch 1.9.0.

	\subsection{Test settings}

	In our numerical experiments, to illustrate the practical performance of \ref{NEPen}, we choose the well-known Fashion-MNIST \cite{xiao2017FashionMinst} and EMNIST \cite{cohen2017emnist} image classification datasets. We develop the network based on the PyTorch examples on image classification. The detailed structure of the tested network is presented in Figure \ref{fig:fcn32}, where
	we enforce the orthogonality of the weights in the first full-connected layer. Let the weight of the first full-connected layer and all the other weights be denoted as $W_1 \in \bb{R}^{1936\times 128}$ and  $W_2$, respectively, then we train the neural network with the constraints $W_1\tp W_1 = I_{128}$. Therefore, training this OCNN can be expressed as solving the following optimization problem,
	\begin{equation}
		\label{Obj_NN}
		\min_{W_1, W_2} ~ f_{NN}(W_1, W_2), \qquad \text{s. t. } W_1\tp W_1 = I_{128}.
	\end{equation}

	\begin{figure}[htbp]
		\centering
		\includegraphics[width=0.8\linewidth]{fcn32}
		\caption{Structure of the OCNN. The neural network has $3$ convolution layers. An average pooling and $3$ fully connected layers follow the convolution layers, where we enforce the orthogonality of the weights in the first fully connected layer.  The negative log-likelihood loss is used as the loss function.}
		\label{fig:fcn32}
	\end{figure}

	Then based on the formulation of \ref{NEPen}, we can reshape \eqref{Obj_NN} as the following unconstrained optimization problem
	\begin{equation}
		\min_{W_1, W_2} ~ f_{NN}(\A(W_1), W_2) + \frac{\beta}{4} \norm{W_1\tp W_1 - I_p}\fs,
	\end{equation}
	which can be solved by existing unconstrained solvers provided in the PyTorch platform.

	Moreover, we compare the proposed approach with existing trivalization approaches \cite{lezcano2019trivializations}. Based on matrix exponential mapping, \cite{lezcano2019trivializations} reformulates \eqref{Obj_NN} as
	\begin{equation}
		\min_{Z_1, W_2} ~ f_{NN}(\mathrm{exp}(Z_1)\tilde{Z}, W_2), \qquad \text{s. t. } Z_1 \in \bb{R}^{1936\times 1936}, Z_1 +  Z_1\tp = 0,
	\end{equation}
	where $\mathrm{exp}$ is the matrix exponential and $\tilde{Z} \in \ca{S}_{1936\times 128}$ is a randomly prefixed matrix to ensure that $\mathrm{exp}(Z_1)\tilde{Z}$ is orthogonal.
	However, as mentioned in \cite{helfrich2018orthogonal, manton2002optimization}, computing the  exponential function for matrices is extremely costly. Therefore, they proposed an improved trivalization approach based on Cayley transform,
	\begin{equation}
		\min_{Z_1, W_2} ~ f_{NN}(\mathrm{cay}(Z_1)\tilde{Z}, W_2), \qquad \text{s. t. } Z_1 \in \bb{R}^{1936\times 1936}, Z_1 +  Z_1\tp = 0,
	\end{equation}
	where $\mathrm{cay}(Z_1) := (I_{1936} + Z_1)^{-1}(I_{1936} - Z_1) $ and $\tilde{Z} \in \ca{S}_{1936\times 128}$ is a prefixed matrix on the Stiefel manifold. Although these approaches require computing the exponential function or inverse of square matrices, they transform \eqref{Obj_NN} into an unconstrained optimization problem
	which can be solved directly by existing unconstrained solvers provided by
	the widely-used PyTorch platform .

	Recently, a collection of PyTorch-based solvers, named McTorch \cite{meghwanshi2018mctorch},  provides pre-defined full-connected and convolution layers. Therefore, we could define the model \eqref{Obj_NN} on PyTorch platform based on McTorch and employ several specialized Riemannian solvers to train it. However, based on the framework provided in \cite{Absil2009optimization}, transforming existing unconstrained solvers to their Riemannian versions requires deep modification to the corresponding unconstrained solvers.   Therefore, only Riemannian stochastic gradient method, Riemannian version of adagrad, and a specialized Riemannian adaptive stochastic gradient method are provided in McTorch. Various other efficient solvers, including Adam \cite{kingma2014adam}, Nesterov accelerated Adam \cite{dozat2016incorporating}, Lamb \cite{you2019large}, Lookahead \cite{zhang2019lookahead}, and AdamW \cite{loshchilov2017decoupled} are still absent for existing Riemannian optimization packages.

	\subsection{Numerical results}
	In our numerical tests, we compare the numerical performance of training orthogonally constrained CNN by \ref{NEPen}, trivialization by exponential mapping (Exp) \cite{lezcano2019trivializations}, trivialization by Cayley transformation (Cayley) \cite{helfrich2018orthogonal, manton2002optimization}, nonsmooth penalty function (L1), quadratic penalty function (Courant) and Riemannian solvers from McTorch package. We choose the penalty parameter $\beta = 0.01$ in \ref{NEPen}, L1 and Courant, then we train the orthogonally constrained CNN by PyTorch build-in unconstrained solvers Adam, AdamW, NAdam, and SGD. Furthermore, we also test the Lamb solver from the torch-optimizer package, which provides various extended solvers for PyTorch.
	However, among all the aforementioned solvers, only Riemannian SGD is provided in McTorch package.
	Therefore, we only use the Riemannian SGD from the McTorch package to train the orthogonally constrained CNN in \eqref{Obj_NN}.
	 For Adam, NAdam, AdamW, and Lamb, we set the learning rate as $0.0005$ while  fix all the other parameters as default. For SGD, we set the learning rate as $0.05$, the momentum factor as $0.9$, and fix all the other parameters as default. Furthermore, we decay the learning rate of each parameter group by $0.8$ after every epoch.

	All the experiments are repeated $5$ times, and the batch size is set as $64$.
	Table \ref{Table_compare} presents the averaged performance of the
	selected approaches in training the proposed network with different datasets and different solvers.
	In Table \ref{Table_compare}, compared with the existing penalty function approaches, \ref{NEPen} achieve similar training speed. Moreover, the solutions obtained by minimizing \ref{NEPen} has much higher accuracy in feasibility than the solutions from other
	penalty-based approaches, while they all have comparable test accuracy.

	Moreover,  compared with the existing trivialization approaches, \ref{NEPen} achieve similar accuracy as the trivialization approaches by an exponential function and Cayley transform. However, as computing the matrix exponential function and matrix inverse are expensive in GPU, those trivialization approaches take longer CPU time than \ref{NEPen} significantly. Moreover, since GPU utilizes  half-precision floating-point numbers, training in GPU by trivialization approaches introduces higher round-off errors, resulting in lower feasibility of the weights in practice.

	Furthermore, compared with existing Riemannian solvers, \ref{NEPen} achieves similar accuracy. However, training orthogonally constrained neural networks by \ref{NEPen} only requires matrix-matrix multiplication, while  the solvers in McTorch compute the retractions in each iteration by singular value decomposition or QR factorization. Therefore, running SGD to minimize \ref{NEPen} can be much faster than the Riemannian SGD provided by McTorch.  Moreover, in the presence of difficulties in developing efficient Riemannian solvers, existing PyTorch-based Riemannian solvers are limited to SGD and Adagrad. In contrast, training orthogonally constrained neural networks by \ref{NEPen} is highly adaptive to existing unconstrained solvers, and various existing efficient solvers can be applied to minimize \ref{Prob_Ori} through \ref{NEPen}.

	\begin{table}[]
		\scriptsize
		\centering
		\begin{tabular}{c@{\extracolsep{\fill}}c|ccc|ccc|ccc@{\extracolsep{\fill}}}
			\hline
			\multicolumn{2}{c}{\multirow{2}{*}{Test instances}} &\multicolumn{3}{@{}c@{}}{Fashion-MNIST} & \multicolumn{3}{@{}c@{}}{EMNIST-Letters} & \multicolumn{3}{@{}c@{}}{EMNIST-Balanced} \\
			\cmidrule{3-5}\cmidrule{6-8} \cmidrule{9-11}
			&  & Acc    & Feas    &  Time (s)/epoch    & Acc    & Feas    & Time (s)/epoch      & Acc    & Feas    &  Time (s)/epoch         \\ \hline
			\multirow{6}{*}{ADAM} & \ref{NEPen}  & 92.79 & 1.68e-04 & 14.85 & 94.47 & 1.37e-04 & 29.99 & 88.73 & 1.60e-04 & 27.13         \\
			& Exp  & 92.75 & 9.03e-03 & 195.63 & 94.60 & 1.37e-02 & 403.23 & 88.71 & 1.43e-02 & 364.84 \\
			& Cayley  & 92.37 & 3.65e-03 & 109.76 & 94.47 & 9.65e-03 & 229.44 & 88.84 & 1.04e-02 & 206.92 \\
			& McTorch &    -    &     -    &     -   &    -     &       -  &       -  &    -      &-         &    -     \\
			& L1  & 92.69 & 2.54e-02 & 14.68 & 94.24 & 5.16e-02 & 30.10 & 88.71 & 5.41e-02 & 27.39       \\
			& Courant  & 92.50 & 1.38e+01 & 14.55 & 94.21 & 2.02e+01 & 30.09 & 88.46 & 2.14e+01 & 27.04      \\ \hline
			\multirow{6}{*}{Adam-W}
			& \ref{NEPen}  & 92.74 & 3.13e-04 & 14.40 & 94.53 & 1.07e-04 & 29.84 & 88.84 & 1.83e-04 & 26.79  \\
			& Exp  & 92.87 & 9.58e-03 & 196.27 & 94.39 & 1.38e-02 & 405.06 & 88.84 & 4.00e-03 & 358.81 \\
			& Cayley  & 92.33 & 3.87e-03 & 110.18 & 94.45 & 1.04e-02 & 228.32 & 88.90 & 1.08e-02 & 209.47 \\
			& McTorch &    -    &     -    &     -   &    -     &       -  &       -  &    -      &-         &    -     \\
			& L1  & 92.51 & 2.01e-02 & 14.78 & 94.08 & 1.53e-02 & 30.17 & 88.37 & 2.48e-02 & 27.20 \\
			& Courant  & 92.08 & 1.29e+01 & 14.48 & 94.15 & 1.89e+01 & 29.76 & 88.21 & 2.05e+01 & 26.83 \\ \hline
			\multirow{6}{*}{SGD}
			& \ref{NEPen}  & 91.50 & 6.26e-05 & 14.40 & 93.94 & 5.57e-05 & 29.50 & 87.89 & 1.14e-04 & 26.64 \\
			& Exp  & 91.46 & 2.30e-03 & 165.16 & 93.26 & 2.20e-03 & 343.20 & 87.49 & 2.27e-03 & 311.08 \\
			& Cayley  & 89.61 & 2.57e-03 & 108.58 & 93.28 & 2.38e-03 & 226.84 & 86.98 & 2.35e-03 & 204.88 \\
			& McTorch   & 91.43 & 2.72e-06 & 23.37 & 93.89 & 2.68e-06 & 48.30 & 87.92 & 2.60e-06 & 44.17 \\
			& L1  & 90.86 & 1.37e-02 & 14.78 & 93.67 & 1.13e-02 & 30.07 & 87.31 & 1.29e-02 & 27.29 \\
			& Courant   & 89.61 & 7.20e+00 & 14.52 & 92.93 & 7.42e+00 & 30.04 & 86.05 & 7.50e+00 & 26.90 \\ \hline
			\multirow{6}{*}{Lamb} & \ref{NEPen}   & 92.96 & 6.87e-05 & 14.43 & 94.45 & 9.16e-05 & 30.07 & 88.65 & 4.58e-05 & 26.99 \\
			& Exp & 93.10 & 2.34e-03 & 170.74 & 94.26 & 2.52e-03 & 354.44 & 88.58 & 2.59e-03 & 322.00 \\
			& Cayley  & 92.42 & 2.61e-03 & 112.94 & 94.40 & 2.46e-03 & 235.13 & 88.63 & 2.56e-03 & 207.35 \\
			& McTorch &    -    &     -    &     -   &    -     &       -  &       -  &    -      &-         &    -     \\
			& L1   & 93.02 & 1.47e-04 & 14.80 & 94.22 & 6.22e-03 & 30.58 & 88.62 & 4.73e-03 & 27.85 \\
			& Courant  & 91.98 & 7.06e+00 & 14.59 & 94.28 & 6.89e+00 & 31.41 & 88.64 & 6.66e+00 & 27.36 \\ \hline
			\multirow{6}{*}{NAdam} & \ref{NEPen}   & 92.49 & 8.39e-05 & 14.46 & 94.52 & 1.98e-04 & 29.68 & 88.86 & 3.05e-04 & 26.80 \\
			& Exp & 92.72 & 2.12e-02 & 191.58 & 94.41 & 4.43e-02 & 403.93 & 88.23 & 4.60e-02 & 356.74 \\
			& Cayley  & 92.70 & 9.59e-03 & 111.35 & 94.45 & 6.88e-03 & 229.76 & 88.89 & 7.06e-03 & 201.21 \\
			& McTorch &    -    &     -    &     -   &    -     &       -  &       -  &    -      &-         &    -     \\
			& L1  & 92.46 & 4.61e-02 & 14.76 & 94.13 & 7.17e-02 & 30.14 & 88.50 & 2.57e-02 & 27.47 \\
			& Courant  & 92.68 & 9.27e+00 & 14.42 & 94.15 & 1.58e+01 & 29.78 & 88.48 & 1.45e+01 & 27.06 \\ \hline
		\end{tabular}
		\caption{A comparison between \ref{NEPen} and other approaches on training orthogonally constrained neural networks. ``Acc'' refers to the test accuracy of the final results. ``Feas'' denotes the feasibility of the final weights in the orthogonally constrained layer. ``Time (s)/epoch'' refers to the averaged running time for each epoch.   }
		\label{Table_compare}
	\end{table}

	\section{Conclusion}

	Nonsmooth nonconvex optimization problems over the Stiefel manifold have real-world applications in various areas.
	The difficulties in handling such problems often arise from  the failure of the chain rule once the required regularity conditions of the objective are missing.
	We present a novel constraint dissolving function \ref{NEPen} for  optimization problems over the Stiefel manifold \ref{Prob_Ori}, whose objective functions are only assumed to be locally Lipschitz continuous. We show that \ref{Prob_Ori} and \ref{NEPen} share the same first-order stationary points and local minimizers in a neighborhood of the Stiefel manifold. The relationships between \ref{Prob_Ori} and \ref{NEPen} illustrate that minimizing a locally Lipschitz smooth function over the Stiefel manifold can be transferred into minimizing \ref{NEPen} without any constraint. Moreover, we show that the Clarke subdifferential of \ref{NEPen} is easy to achieve from $\partial f$. Therefore, the exactness and accessibility of Clarke subdifferential enable the direct implementation of various unconstrained optimization approaches for solving \ref{Prob_Ori} through \ref{NEPen}.

	We present a representative example to demonstrate that \ref{NEPen} admits direct implementation of existing unconstrained optimization approaches, while their theoretical properties are simultaneously retained.
	We then introduce a generalized framework for a class of subgradient methods and establish their convergence theories based on \cite{davis2020stochastic}.
	Furthermore, numerical examples on an orthogonally constrained
	neural network for classification problems
	demonstrate that the problem can solved via \ref{NEPen} by efficient and direct implementation of existing unconstrained  optimization solvers.