Find $\mathbf{A}$ minimizing $\lVert \mathbf{A}^{\top}\mathbf{AX}-\mathbf{X} \rVert _{F}^{2}$ with a given $\mathbf{X}$ $1.$ Problem:
There are three subproblems:
$(1.1)$ Given a data matrix $\mathbf{X}\in \mathbb{R}^{N\times D}$, find a matrix $\mathbf{A}\in \mathbb{R}^{M\times N} (M\ll N)$ satisfying:
$$
\mathbf{A}=\underset{\mathbf{A}}{\arg\min}\lVert\mathbf{A}^\top\mathbf{A}\mathbf{X}-\mathbf{X}\rVert_F^2.
$$
$(1.2)$ Given a data matrix $\mathbf{X}\in \mathbb{R}^{N\times D}$, find a bipolar matrix $\mathbf{B}\in \mathbb{Z}^{M\times N} (M\ll N)$ and a scaling factor $\alpha\in\mathbb{R}$ satisfying:
$$
\mathbf{C}=\underset{\mathbf{C}}{\arg\min}\lVert\mathbf{C}^\top\mathbf{C}\mathbf{X}-\mathbf{X}\rVert_F^2,~~~~\text{s.t.}~\mathbf{C}=\alpha\mathbf{B},
$$
$~~~~~~$where each element of $\mathbf{B}$ can only be choosed from $\{-1,+1\}$.
$(1.3)$ Based on $(1.2)$, could you please get an optimal binary matrix $\mathbf{C}$ (i.e., each element of $\mathbf{B}$ can only be choosed from $\{0,1\}$)?

$2.$ Background:
(Compressed Sensing)
$(2.1)$ I project an $N$-dimensional dataset (with $D$ data points) into an $M$-dimensional space by $\mathbf{Y}=\mathbf{AX}$ and then reconstruct it by $\mathbf{\hat{X}}=\mathbf{A^\top Y}=\mathbf{A^\top AX}$. I want to get an optimal $\mathbf{A}$ such that the $\ell_2$-error between $\mathbf{X}$ and $\mathbf{\hat{X}}$ is minimized.
$(2.2)$ For higher compression efficiency and lower memory complexity of $\mathbf{A}$, I want to get a bipolar $\mathbf{C}$ that may works well or close to floating-point $\mathbf{A}$.

$3.$ My Efforts:
$(3.1)$ I think I can solve subproblem $(1.1)$ by singular value decomposition (SVD) of $\mathbf{X}=\mathbf{U\Sigma V^\top}$ with singular values $\sigma_i=\mathbf{\Sigma}_{i,i}$ satisfying $\sigma_1\ge \sigma_2\ge\cdots\ge \sigma_N$. The optimal $\mathbf{A}$ should be formed by the first $M$ rows of $\mathbf{U}^\top$ [Reference].
$(3.2)$ For subproblem $(1.2)$, I guess that the optimal scaling factor $\alpha=1/\sqrt{N}$ since the optimal $\mathbf{C}$ should be row-normalized, and I do not know how to get the optimal $\mathbf{B}$.
$(3.3)$ I guess that I can get a bipolar $\mathbf{B}$ by discretizing each element of $\mathbf{A}$ with
$$
b_{i,j}=f\left( a_{i,j} \right) =\left\{ \begin{array}{l}
 -1,\ a_{i,j}<0\\
 1,\ a_{i,j}\ge 0\\
\end{array} \right.,
$$
$~~~~~~~~~$which may be better than random ones (by Bernoulli sampling) but not optimal?

$4.$ Some Other Minor Questions:
$(4.1)$ I do not know if I can directly write that “$\mathbf{B}\in\{-1,+1\}^{M\times N}$”?
 A: $
\def\b{\beta}\def\l{\lambda}\def\s{\sigma}
\def\o{{\tt1}}\def\p{\partial}
\def\LR#1{\left(#1\right)}
\def\size#1{\operatorname{size}\LR{#1}}
\def\trace#1{\operatorname{Tr}\LR{#1}}
\def\qiq{\quad\implies\quad}
\def\grad#1#2{\frac{\p #1}{\p #2}}
\def\gradd#1#2{\frac{d #1}{d #2}}
\def\c#1{\color{red}{#1}}
\def\CLR#1{\c{\LR{#1}}}
\def\gradLR#1#2{\LR{\grad{#1}{#2}}}
\def\fracLR#1#2{\LR{\frac{#1}{#2}}}
$You can use the logistic function to find an approximate solution to subproblem $(1.3)$
$$\eqalign{
\s(z) &= \fracLR{\o}{\o+e^{-z}} &\qiq \gradd{\s}{z} = \s-\s^2  \\
}$$
Introduce an unconstrained matrix $U$ as the new independent variable and define $C$ as the element-wise logistic function
$$\eqalign{
C &= \s(U) = {J}\oslash\LR{J+e^{-U}} \qiq\c{dC = (C-C\odot C)\odot dU} \\
}$$
where $J$ is the all-ones matrix and $(\odot,\oslash)$ denote  elementwise multiplication and division.
Use the gradient wrt $C$ to compute the gradient wrt $U$
$$\eqalign{
\phi &= \|C^TCX-X\|_F^2 \\
\grad{\phi}{C} &= 2CC^TCXX^T + 2CXX^TC^TC - 4CXX^T \;\doteq\; H\\
d\phi &= H:\c{dC} \\
  &= H:\CLR{(C-C\odot C)\odot dU} \\
  &= \LR{(C-C\odot C)\odot H}:dU \\
\grad{\phi}{U} &= (C-C\odot C)\odot H \;\doteq\; G \\
}$$
where $(:)$ denotes the Frobenius product, which is a concise
notation for the trace
$$\eqalign{
A:B &= \sum_{i=1}^m\sum_{j=1}^n A_{ij}B_{ij}
 \;=\; \trace{A^TB} \\
A:A &= \|A\|^2_F \\
A:(B\odot C) &= \sum_{i=1}^m\sum_{j=1}^n A_{ij}B_{ij}C_{ij}
 \;=\; (A\odot B):C \\
\\
}$$

Now that you have an analytical expression for the gradient, you can use
any gradient descent method to solve the zero gradient condition.
A typical step in such an algorithm will look like
$$\eqalign{
U_{k+1} &= U_k - \l_k G_k \\
}$$
where $k$ is the iteration counter and $\l_k$ is the step-length parameter which will depend on the particular gradient method that is used. Personally, I find the Barzilai-Borwein method using Raydan's non-monotone line search to be effective and easy to code.
Using the vanilla logistic function, the components of $A$ are approximately equal to zero or one, but this approximation can be improved by introducing a sharpening parameter $\b\gg\o$
$$\eqalign{
\s_\b(z)
 &= \fracLR{\o}{\o+e^{-\b z}}
 &\qiq \gradd{\s_\b}{z} = \b\LR{\s-\s^2}  \\
}$$
As the parameter increases, the components of $A$ converge to either zero or one.
