I believe that there is a good alternative that benefits from the geometry of the parameter space and completely eliminates the need for constrained optimization. If you explicitly wanted to make use of Lagrangians, I will definitely not be answering the question, but I thought it might be worthwhile to consider the perspective I will describe.
In particular, the approach I will be presenting uses Riemannian manifold optimization tools.
Note that the parameter of interest $\mathbf{X}\in \mathbb{R}^{n \times k}$ is simply an element of the Stiefel Manifold $\mathcal{M}\equiv V_k(\mathbb{R}^n)$ (the set of all orthonormal $k$-frames in $\mathbb{R}^n$):
$$
\begin{align}
\mathbf{X} &\in V_k(\mathbb{R}^n) \subset \mathbb{R}^{n\times k}\\
V_k(\mathbb{R}^n) = \{\mathbf{X}&\in\mathbb{R}^{n\times k} : \mathbf{X}^\top\mathbf{X}=\mathbf{I}\}.
\end{align}
$$
With this definition, we can benefit from the Riemannian structure of $V_k(\mathbb{R}^n)$ and optimize the energy with any unconstrained optimization algorithm, such as the Riemannian-gradient descent:
\begin{align}
&\text{ While } \mathbf{X}_{k} \text{ does not sufficiently minimize } f \\
&\text{$\quad$- Pick a gradient related descent direction }\boldsymbol{\eta}_k\in T_{\mathbf{X}_k}\mathcal{M}\\
&\text{$\quad$- Choose a retraction } R_{\mathbf{X}_k}:T_{\mathbf{X}_k}\mathcal{M}\rightarrow \mathcal{M}.\\
&\text{$\quad$- Choose a step length } \tau_k\in \mathbb{R}.\\
&\text{$\quad$- Set } \mathbf{X}_{k+1}\gets R_{\mathbf{X}_k}(\tau_k\boldsymbol{\eta}_k).\\
&\text{$\quad$- } k\gets k+1
\end{align}
Here $f$ denotes the function we optimize:
$$
f = \frac{1}{2}\|\mathbf{X}\mathbf{X}^\top - \mathbf{A} \|_\mathcal{F}^2 = \frac{1}{2}\text{tr}\big((\mathbf{X}\mathbf{X}^\top-\mathbf{A})(\mathbf{X}\mathbf{X}^\top-\mathbf{A})^\top\big)
$$
whose gradient reads:
$$
\nabla_{\mathbf{X}}f = (\mathbf{X} \mathbf{X}^\top -\mathbf{A})\mathbf{X}+(\mathbf{X}\mathbf{X}^\top +(-\mathbf{A})^\top )\mathbf{X}.
$$
We generally pick the gradient related direction $\boldsymbol{\eta}$ as the projection of the Euclidean gradient onto the tangent space of the manifold (or for a broad class of manifolds including the Stiefel manifold, we can instead use the logarithmic-map):
$$
\boldsymbol{\eta}\triangleq \text{grad }_{\mathbf{X}} f = \Pi_{\mathbf{X}} \Big( -\nabla_{\mathbf{X}} f \Big)
$$
For now let us assume that $\tau_k$ is fixed, e.g. $\tau_k=0.1$.
The retraction operator $R(\cdot)$ of the Stiefel manifold is analytically defined or in other words the true exponential map is available. More details on that are given in this math-se post and of course in the seminal paper of Edelman et. al:
Edelman, Alan, Tomás A. Arias, and Steven T. Smith. "The geometry of
algorithms with orthogonality constraints." SIAM journal on Matrix
Analysis and Applications 20.2 (1998): 303-353.
All the Riemannian operators for Stiefel manifold are included in the toolboxes such as Manopt or ROPTLIB. Finally, we can update the obtain the new value $\mathbf{X}_{k+1}$ by:
\begin{align}
\mathbf{X}_{k+1} = R_{\mathbf{X}_k}(\tau_k \boldsymbol{\eta}_k)
\end{align}
Note that $\mathbf{X}_{k+1}$ will be an orthonormal $k$-frame obeying the aforementioned constrained naturally.
Riemennian gradient descent is the simplest choice of Riemannian optimization and there are many others such as Riemannian-LBFGS or Riemannian-Trust Region. Many of those choices vary in how they compute the step size $\tau_k$, usually by a form of line-search.
More recently, Hu et al. considered a very similar problem to the one of this questions where the minimization with orthonormality constraints (again, Stiefel manifold) is made efficient. This method too, uses the Riemannian structure of the problem.
