We observe that NMF-LR achieves the smallest reconstruction error among all methods but suffers for the classification task.  This is expected from the construction of the dataset, as the synthetic images are nonnegative linear combinations of images of digits '2' and '5', and the same dictionary atoms are also used as filters for label prediction.  On the other hand, logistic regression on pixels does not compress the data matrix so we assigned a relative reconstruction error of 1. In Figure \ref{fig:MNIST_pareto}, we observe that, except in a few instances, most SDL models lie between the two extremes of (1) LR and (2) NMF-LR in the sense that achieving significantly better classification accuracies (both in terms of accuracy and F-score) with small reconstruction error. For instance, SDL-filt with $ξ=10$ achieves a relative reconstruction error of about 0.22 and classification accuracy above 80$%$, which is more than double the performance of NMF-LR and also about $11%$ better than LR accuracy. It is interesting to note that SDL-conv-filt with $ξ=0.1$ achieves the best classification accuracy of about $91%$ but its reconstruction error is quite large at around $0.7$.

For more qualitative analysis, we plot various estimated dictionaries $$  and compare them. Figure \ref{fig:MNIST_filt_dict} shows how the dictionary matrix $$ estimated by SDL-filter changes depending on the level of tuning parameter $ξ∈{0.01, 0.1, 1}$. When $ξ=1$, the combined SDL loss in \eqref{eq:ASDL_1} puts some significant weight on the matrix factorization term, so the learned dictionary $$ should be reconstructive of the synthesized images in $_$. Indeed, the learned atoms in Figure  \ref{fig:MNIST_filt_dict} left shows shapes of digits `5' and `2'. Further, the second atom resembling '2' is associated with a negative regression coefficient, indicating that being close to `2' may be partially aligned with being close to `7', which corresponds to being a negative example. On the other hand, decreasing the value of $ξ$ increases the amount of supervision. The learned atoms in Figure \ref{fig:MNIST_filt_dict} middle and right resembles less of the digits '2' and '5', but it seems that some abstract shape with large positive (for $ξ=1$) and large negative (for $ξ=0.01$) are learned and the resulting classification accuracies increase. The other atoms seem to be the shape of `8', which are seemingly learned from superpositions` 2' and `5'. One can regard the learned dictionary atoms as the `best effort' of SDL-filter in balancing the two partially aligned reconstructions and discrimination taken from the space of images of linear combinations of `2' and `5'.

%These three $\hat{W}_x$ were obtained by the filter-based method. When $\xi=0$, the loss function does not contain matrix factorization loss, therefore, the $\hat{W}_x$ has highlighted pixels rather than the shape of a specific number related to the class (left). On the contrary, when $\xi=1$, $\hat{W}_x$ is very close to the image of original basis matrix $W_{{true}_X}$ (right). $\hat{W}_x$ with $\xi = 0.5$ was obtained as the combination of logistic regression loss and matrix factorization loss with equal weight (middle). The second basis matrix looks 8 generated by the combination of 2 and 5 properly to predict the second class, 8.

\begin{figure}[h!]
\centering
\includegraphics[width=1\linewidth]{Figures/MNIST_filt_dict_xi.pdf}
\caption{Estimated dictionary matrix $\hat{\W}$ from SDL-filter  depending on the level of tuning parameter $\xi\in \{ 1,0.1,0.01 \}$. The smaller the tuning parameter is the stronger the supervision effect is.}
\label{fig:MNIST_filt_dict}
\end{figure}

\begin{figure}[h!]
\centering
\includegraphics[width=1\linewidth]{Figures/MNIST_dict_comparison.pdf}
\caption{Estimated basis matrix $\hat{W}_x$ from filter based method depending on the level of tuning parameter $\xi$. $\xi = 0$ (Left), $\xi = 0.5$ (Middle), $\xi = 1$ (Right)}
\label{fig:MNIST_dict_comparison}
\end{figure}

In Figure \ref{fig:MNIST_dict_comparison}, we plot the estimated  dictionaries $$ from all five methods except LR with chosen level of tuning parameter $ξ$. It is interesting that NMF learned almost identical atoms (with inferior classification performance) that resemble the typical shape of nonnegative linear combinations of digits `2' and `5', instead of learning separate atoms of shapes `2' and `5' individually.  For convex SDL algorithms, we find that typically a larger tuning parameter is required for fast convergence, which we also observed in Figure \ref{fig:benchmark_MNIST}.

\section{Applications}
\subsection{Supervised Topic Modeling on fake job postings dataset}

%Fake job postings have been causing financial losses to individuals seeking jobs.

According to Better Business Bureau, a non-profit organization that monitors and evaluates job postings, there were 3,434 fake job postings reported in 2019. These scam postings result in huge financial loss and the average loss per victim is \$3,000 according to the FBI reports [7]. In this section, we use our methods to classify fake job postings on a dataset 'real-or-fake-job posting-prediction' in Kaggle [6]. In doing so, the new method can also simultaneously learn topics that are most effective in classifying fake job postings.	We compared the performance of classifiers based on multiple indexes such as accuracy, and F-score, and find what factors are highly associated with fake job postings. Identifying the relevant characteristics of scams and building prediction models will help prevent potential financial losses in advance.

\subsubsection{Dataset description and preliminary analysis}
\label{sec:fakejob_setup}

There are 17,880 postings and 15 variables in the dataset including binary variables, categorical variables, and textual information of \textit{job description}. Among the 17,880 postings, 17,014 are true job postings (95.1\%) and 866 are fraudulent postings (4.84\%), which shows a high imbalance between the two classes. This imbalance is the main characteristic of the dataset. We coded fake job postings as positive examples and true job postings as negative examples. Due to the high imbalance, the accuracy of classification can be trivially high (e.g., by classifying everything to be negative), and hence achieving a high F-score is of importance.

In our experiments, we represented each job posting as a $p=2480$ dimensional word frequency vector computed from its \textit{job description} and augmented with $q=72$ auxiliary covariates of binary and categorical variables including indicators of the posting having a company logo or the posted job  in the United States or not. For computing the word frequency vectors, we represent the job description variable as a term/document frequency matrix with Term Frequency-Inverse Document Frequency (TF-IDF) normalization [60]. TF-IDF measures the relative importance of each term in a collection of documents. If a word is common to all documents, then it is less likely to have an important meaning. %After filtering,
The top 2480 most frequent words were used for the analysis. %The below Figure \ref{fig:Figure4} shows the 10 most frequent words from the variable according to fraudulent and true job postings.

\begin{figure}[h!]
\centering
\includegraphics[width=1\textwidth]{Figures/fakejob_LR_coefficients.pdf}
\caption{The top 20 variables with largest coefficients from logistic regression on $p=2480$ words in job description (left) and $p+q=2480+72$ words and auxiliary covariates combined (right). In the left panel, red bars indicate the words that appear as dominant keywords in the forthcoming topic modeling analysis. In the right panel, blue bars indicate words and grey bars indicate auxiliary covariates.}
\label{fig:fakejob_LR}
\end{figure}

As a preliminary analysis, we first apply logistic regression either on the $p$-dimensional word frequency vectors or on the $(p+q)$-dimensional combined feature vectors (Figure \ref{fig:fakejob_LR} right). For the former experiment, Figure \ref{fig:fakejob_LR} left shows 10 words each with positive and negative regression coefficients with the largest absolute values. The results indicate that having a large frequency of words such as `earn', `money', and `link' is positively correlated with being a fake job, whereas postings with a high frequency of words such as `team', `client', and `website' as well as with company logo and jobs outside of the US are more likely to be true job postings.

\subsubsection{Supervised Topic Modeling with Auxiliary Covariates}

Topic modeling is a classical technique in text data analysis that seeks to find a small number of `topics', which are groups of words that share semantic context. The grounding assumption is that a given text may be built upon such topics as latent variables. Methods such as nonnegative matrix factorization (NMF) [40] and latent Dirichlet allocation (LDA) [65, 15, 33] have been successfully used to detect or estimate such latent semantic factors. Also, `supervised' topic modeling techniques have been studied, where one seeks to group words not only by their semantic contexts but also using their `functional contexts' that are provided by additional class labels. See, for example, [47] for LDA-based approaches and [29] for NMF combined with linear regression model (see \eqref{eq:SDL_regression_H1}).
Here we mainly compare two methods, namely, (1) NMF with logistic regression and (2) SDL-filter with nonnegativity constraints on $$ and $$̋. However, we do compare the performance of all four SDL models in Figures <ref> and <ref>. We note that for the purpose of topic modeling, it is crucial to use nonnegativity constraint on the dictionary matrix $\W$ in the SDL model (<ref>) as word frequencies are nonnegative and we would like to decompose a given document's word frequency as additively rather than subtractively in order for better interpretability (see, e.g., [40]).

Comparison between (supervised) topics learned by NMF and SDL-filter for the fake job posting data. (a) Four out of 25 topics learned by NMF are shown together with the corresponding logistic regression coefficients. (b) Four out of 20 supervised topics learned by SDL-filter with tuning parameter $\xi=5$ are shown together with the corresponding logistic regression coefficients. (c) Similar as (b) but with tuning parameter $\xi=1$. (d) Nine out of 20 supervised topics (white background)+ 72 auxiliary covariates (dark background) learned from SDL-filt with 72 auxiliary covariates are shown with their corresponding logistic regression coefficients. Corresponding classification accuracy and F-scores are also shown in the subtitles with fake job postings being the positive examples. For topic wordclouds (white background), word sizes are proportional to their frequency.

First consider  Figure <ref> (a), which shows topics (shown as wordclouds) learned by NMF and their associated regression coefficients. Namely, after learning a dictionary matrix $\W\in\R^{p\times 25}$ by NMF from the job description matrix of shape $p\times n$ with $n=17,780$, each of the $r$ columns of $\W$ becomes the topic frequency vector and top 10 words with higheset frequencies are shown as wordcloud. NMF was able to find topics that summarize specific job information. More specifically, the upper right and lower left topics correspond to beauty and healthcare-related jobs. However, as can be seen by the low F-score reported in Figure <ref> (a), while the 25 topics learned by NMF give generic job descriptions, they may not be helpful to determine if a job posting is indeed fake. The main reason that we have these topics is that the dataset is highly imbalanced. Since most of the postings are true job postings (95%), when we first conduct dimension reduction based on NMF, the topics that we learned are mainly determined by dominant true job postings, rather than fake job postings.

On the other hand, some selected supervised topics (out of $20$ total) learned by SDL-filter with $\xi=5$ and $1$ are shown in Figures <ref> (b) and (c). In each case, the upper left and lower right topics are the ones with the largest positive and negative regression coefficients, respectively, and the upper right and lower left ones are manually selected for illustration purposes. For $\xi=1$ in Figure <ref> (c), notice that the upper left topic with positive regression coefficient consists of words that appear frequently on fake job postings (e.g., `money', `earn', `link'), while the lower right topic uses the words from true job postings (e.g., `team', `client', `website'), both detected by logistic regression in Figure <ref>. Topics with neutral regression coefficients are mainly used to reconstruct data matrix rather than for classification purposes. Note that the corresponding F-score of 0.43 is achieved using only  20 variables (topics) and is on par with the F-scores obtained by logistic regression using $p=2480$ or $p+q=2552$ variables in Figure <ref>.

Increasing the tuning parameter $\xi$ from 1 to $5$ weakens the supervision effect. Accordingly, the two neutral topics in  Figure  <ref> (b) becomes generic job descriptions as found by NMF in Figure  <ref> (a), but the two extreme ones (upper left and lower right) maintain similar content and large absolute values of their regression coefficients.

We also conduct a similar analysis using SDL-filter with $\xi=0.001$ and $r=20$ topics along with 72 auxiliary covariates after converting categorical variables to one-hot-encoding. In Figure <ref> (d), we show the covariates with the largest absolute regression coefficients, which is a mix of supervised topics (white background) and auxiliary variables (dark background). This setting achieves the best F-score of $0.52$ by using $20+72$ variables, which is still significantly less than using all 2552 variables while enjoying better interpretability. We see that SDL-filter automatically combines words that are positively or negatively associated with fake job postings in an ensemble with auxiliary covariates. In other words, SDL-filter seems to perform simultaneous supervised topic modeling on text data, while incorporating auxiliary covariates for improved performance.

§.§.§ Evaluation of model performance

We provide a summary of classification accuracy and F-score of various settings in Table  <ref> (see Tables <ref> and <ref> for the full results). Note that due to the high imbalance in the dataset (only 5% of fake job postings), getting a high classification accuracy is trivial (e.g., by classifying all as true job postings), so getting high $F$-score is more important. One can see that SDL-filter is overall the method of best classification performance, both in terms of accuracy and F-score, which improves when using auxiliary covariates. In contrast to comparable performance of convex SDL algorithms in the semi-synthetic MNIST data in Figure <ref>, for fake job postings dataset they show mediocre  performance. It seems that for larger datasets in high dimension, one needs more extensive hyperparameter tuning for convex SDL methods than the nonconvex ones.

Tables of best average F-score over five runs from each of the six methods for the fake job postings data in Section <ref> with (left) and without (right) the auxiliary covariates for tuning parameter $\xi\in \{0.01, 0.1, 1, 5, 10\}$. See Tables <ref>  and <ref> in Appendix for more details.

Pareto plot of relative reconstruction error vs. classification accuracy/F-score for various models on fake job postings dataset.

As in Figure <ref> for the semi-synthetic MNIST dataset, we also provide a Pareto plot in Figure <ref> to evaluate the performance of various SDL models on the fake job posting dataset against the benchmark models of logistic regression (LR) and NMF followed by logistic regression (NMF-LR). Recall that the Pareto plot shows how a model simultaneously performs two objectives of reducing the reconstruction error $\lVert \X_{\textup{data}} - \W\H \rVert$ as well as increasing the classification accuracy. As before, increasing the tuning parameter $\xi$ in various SDL models seems to interpolate between two extremes of LR and NMF-LR. We observe that SDL-filter performs best overall, in some cases achieving both goals with better classification performance than LR. The inferior performance of convex SDL models on the real dataset in contrast to their superior performance on the semi-synthetic dataset in Figure <ref> indicates that, in practice, the convex SDL models require more hyperparameter tuning. For instance, we did not try to fine-tune the stepsize, which we fixed at $\tau=0.01$ throughout the experiments.

We also report that the topics learned from SDL-feature share similar characteristics with respect to the tuning parameter $\xi$ in figure <ref>. We also mention that the topics learned from convex SDL algorithms, which cannot be used with a nonnegative constraint on the dictionary matrix $\W$, turns out to be uninformative and not very much interpretable. This is expected since the use of nonnegativity constraints on $\W$ and $\H$ (i.e., strong constraints on the SDL model, see \ref{assumption:A1}) is crucial for matrix-factorization-based topic modeling experiments (see [40]). We omit figures for the topics of SDL-feature and the convex SDL models.

%and filter-convex-based dictionary learning have very similar characteristics. From this observation, topics from supervised dictionary learning conduct dimension reduction and classification at the same time, which is different from the usual dimension reduction technique as shown in figure \ref{fig:Figure6} (a). However, topics from feature-convex-based dictionary learning are a collection of unrelated words, neither from general job postings nor from fraudulent job postings. It is mainly because it is based on singular value decomposition for dimension reduction.

%Also, topic modeling provides dimension reduction of the original high-dimensional word frequency representation of text data. In our context, 	instead of using the entire $p=2480$ words, we may learn $r$ (e.g., $r\in \{20,25\}$) topics and work with the reduced $r$-dimensional feature representation for more efficient and effective data analysis. This would also be useful if we want to combine with low-dimensional auxiliary covariates.

% In our context, one may try to find $r$ topics in the job posting dataset and then use the $r$-dimensional compressed representation of the job postings (instead of $p$-dimension) for logistic regression. In this way, we may learn the association between topics and fake/true job postings.

%One way to find \textit{topic} is to use the classic non-negative matrix factorization technique. In other words, first, use non-negative matrix factorization for dimension reduction, and then apply logistic regression using a selected set of words, named topic.

%\yaocomment{It seems that there is no explanation about Table 2. Is that for fake job posting, chest X-ray images, or simulation data?}

\subsection{Supervised dictionary learning on chest X-ray images for pneumonia detection}

Pneumonia is an acute respiratory infection that affects the lungs. According to WHO reports, it accounts for 15\% of all deaths of children under 5 years old, killing 808,694 children in 2017. Moreover, currently, about 15\% of COVID-19 patients suffer from severe pneumonia [73].   Chest X-ray is one inexpensive way to diagnose pneumonia, but rapid radiological
interpretation is not always available. A successful statistical model to classify pneumonia from chest X-ray images will enable rapid pneumonia diagnosis with high accuracy, which will be able to reduce the burden on clinicians and help their decision-making process. In this section, we apply our SDL methods for chest X-ray images for pneumonia detection.

\begin{figure}[h!]
\centering
\includegraphics[width=0.8\textwidth]{Figures/pneumonia_cell.jpeg}
\caption{ The normal chest X-ray (left panel) depicts clear lungs without any areas of abnormal opacification in the image. Bacterial pneumonia (middle) typically exhibits a focal lobar consolidation, in this case in the right upper lobe (white arrows), whereas viral pneumonia (right) manifests with a more diffuse “interstitial” pattern in both lungs. The figure and the description are excerpted from [35].
\label{fig:pneumonia_cell}
\end{figure}

The pneumonia data set was first introduced in [35]. There is a total of 5,863 chest X-ray images from children, consisting of 4,273  pneumonia patients and 1,583 healthy subjects. The images were collected from pediatric patients one to five years old from Guangzhou Women and Children’s Medical Center, Guangzhou. For the analysis of chest X-ray images, all images were initially screened for quality control and two expert physicians diagnosed the images. In the reference [35], an extremely accurate image classification system has been developed (with a classification accuracy of 92.8\%) using sophisticated deep neural network image classifiers. We intend to demonstrate that our SDL methods yield interesting and promising results for medical image classification tasks while being significantly simpler and easier to train than the deep neural network models.

In order to apply our SDL methods, we resize each chest X-ray image into an $180 ×180$ pixel image. Vectorizing each image, we obtain the data matrix $_∈^32,400×5,863$. We label pneumonia images with 1 and normal images with 0, obtaining the labeled matrix $_ ∈{0,1}^1×5,863$. We used the deterministic test/train split provided by the original references [35], where the train and the test sets consist of 5,216 and 624 images, respectively. Standard logistic regression with $p=180^2=32,400$ individual pixels as the explanatory variable yields a classification accuracy of 82\%. However, it is not entirely reasonable to assume individual pixels in the image to be correlated with pneumonia. Instead, we may associate certain latent shapes with pneumonia by using dictionary learning methods.

\begin{figure}[h!]
\centering
\includegraphics[width=1\textwidth]{Figures/pneumonia_dict.pdf}
\caption{25 dictionary atoms learned from chest X-ray images by NMF and SDL-feature with $\xi=50$. Corresponding logistic regression coefficients, as well as classification performances, are also shown. For SDL-feature, we used $L_{1}$ regularization coefficient of $5$ for the code matrix $\H$ and no $L_{2}$-regularization. Positive regression coefficients indicate a positive correlation with having pneumonia. There is a clear contrast between the two extreme atoms (upper left and lower right) according to their correlation with pneumonia.
\label{fig:pneumonia_dict}
\end{figure}

Figure \ref{fig:pneumonia_dict} shows 25 dictionary atoms of size $180\times 180$ learned by NMF (left) and SDL-feature (right) together with their corresponding logistic regression coefficient shown. For SDL-feature, we used tuning parameter $\xi=50$ and $L_{1}$-regularization coefficient of $5$ on the code matrix $\H$. Namely, we used Algorithm \ref{algorithm:SDL} for the feature-based SDL model in \eqref{eq:ASDL_1} together with an additional $L_{1}$ regularization term of $\lambda \lVert \H  \rVert_{1}$ added to the loss function in \eqref{eq:ASDL_1} with $\lambda=5$. Using the $L_{1}$-regularization on the code matrix is standard in dictionary learning literature [50, 52, 49] in order to learn sparse representation on an over-complete dictionary. For this particular application, we find such regularization helps improving the classification accuracy better than using $L_{2}$-regularization on any of the factors. We found the value of $\xi=50$ and $\lambda=5$ by using a grid search for $\xi=\{0.01, 0.1, 1, 10, 50, 100\}$ and $\lambda=\{0,1,5,10\}$.

%\yaocomment{This is for my understanding. Did we add $L_{1}$ regularization in any other numerical studies for SDL models or just for this chest X-ray data? }

Note that the code matrix $\H$ is constrained to be nonnegative during optimization.  Since the probability of the existence of pneumonia is calculated by the logit transformation of $\Beta (\mathbf{1},\H^{T})^{T}$, the positive value of $\Beta$ indicates that the corresponding atom is more related to pneumonia with the magnitude of the regression coefficient indicates the strength of such association, and the negative value of $\Beta$ suggests that the atom is related to normal. Comparing the NMF and SDL-feature dictionaries in Figure \ref{fig:pneumonia_dict}, we find that the regression coefficients for NMF atoms are quite neutral, but there are three atoms (two upper left and one lower right) with an order of magnitude larger absolute regression coefficient learned by SDL-feature. A closer investigation shows that the atom with regression coefficients $-0.066$ has almost no signal (dark) around the lung, whereas the two atoms with regression coefficients 0.052 and 0.046 show the opposite contrast, having most signals (bright) around the lung and weak signals elsewhere. Although one should pay extra caution when interpreting machine learning results as a clinical statement, this observation seems to be well-aligned with some basic characteristics of normal or pneumonia chest X-ray images as shown in Figure \ref{fig:pneumonia_cell} (see also the caption).

Lastly, we report some additional details of these experiments. First, the training of SDL-feature on the chest X-ray dataset and making predictions is extremely efficient and the entire process of training and testing takes under 2 minutes on an average laptop computer (100 iterations on an Apple M1 chip). Second, we find that SDL-filter achieves higher classification performance with accuracy $84.3\%$ and F-score $87.9\%$ with tuning parameter $\xi=0.01$, $r=20$ atoms, and the same $L_{1}$-regularization coefficient on $\H$ of $5$. However, the learned atoms are too noisy and not quite interpretable as the ones learned by SDL-feature in Figure \ref{fig:pneumonia_dict} right. Training and testing take about ten minutes on the same machine. Third, convex SDL models, in this case, take a long time (about an hour) with the same number of iterations and we omit the results.

%with positive regression coefficients and one atom with a negative co

%We conducted the logistic regression with 32,400 explanatory variables, the logistic regression with L1 penalty ($\lambda=10$) and L2 penalty ($\lambda=0.01$), and LMF with 25 latent factors, supervision parameter $\xi=40$ with sparsity constraint at $H$ ($\lambda_{H}=5$). Based on the given train, validation, and test set in Kaggle, we found the optimal $\lambda$ for LASSO and Ridge and evaluated the model performance. The performance of the four models was compared to the result from Convolutional neural network model in the Cell paper. The results given in Table \eqref{tab:Table2} indicate that LMF is comparable with other classical methods.

%In figure \ref{fig:Figure5}, we highlighted regions in yellow that are found to be positively associated with pneumonia based on logistic regression and LMF. More precisely, for logistic regression, the regression coefficient for each pixel in $180\times 180$ image is learned, and the highlighted region consists of pixels with their regression coefficients at least as large as $0.9$ times the maximum regression coefficient. On the other hand, the highlighted regions for LMF consist of pixels at least as large as $0.9$ times the maximum `score', where the score for each pixel is computed by the value in the reconstruction with dictionary images weighted by their regression coefficients. Since logistic regression recognizes a specific region significant, if many patients in the training set have pneumonia in the same particular region, the highlighted region does not vary much regardless of the actual location of pneumonia in each patient. On the other hand, LMF highlights different regions for each patient since its regression coefficient is associated with a geometric pattern, rather than individual pixels.

%{\color{red}This section should contain a clear and convincing discussion of the theoretical properties of the method. It should also discuss under what assumptions the methods should work and under what conditions they will fail.}

\section{Concluding Remarks}

%{\color{red}What is the answer to the question? What did you learn about the methods?}

In this paper, we provided a comprehensive treatment of a large class of supervised dictionary learning methods in terms of model construction, optimization algorithms and their convergence properties, and statistical estimation guarantees for the corresponding generative models. SDL models find best balance between two objectives of data modeling by latent factors and label classification while achieving simultaneous dimension reduction and classification, which makes them suitable for coping with high-dimensional data. We demonstrated that our methods show comparable performance with classical models for classifying fake job postings as well as chest X-ray images for pneumonia while learning basis topics or images that are directly associated with fake jobs or pneumonia, respectively. In addition, the new methods achieve such comparable performance with much reduced number of variables and with much more homogeneous and interpretable models.

Our method can potentially be used for a number of high-dimensional classification problems, especially for areas where interpretability is required such as natural language processing and biomedical image processing. While a number of sophisticated deep-learning-based approaches are gaining popularity due to their extreme success in diverse problems including image classification and voice recognition, an inherent downside is the loss of interpretability due to a severe over-parameterization and sophisticated design of such algorithms. In this work, we showed that combining two classical methods of nonnegative matrix factorization and logistic regression could achieve comparable performance while maintaining the transparency of the method and the interpretability of the results. % and theoretical convergence guarantee.

% AnotherAs another application, the method is applicable to highly unbalanced data, especially for the analysis of rare diseases. Unlike other existing dimensional reduction methods which mainly find latent factors from the normal group, our method jointly uses label information to find a classifier, so we can construct a better-performed model.

One of the main techniques we developed in this work is a `double lifting procedure', which transforms the SDL problem into a CAFE problem \eqref{eq:CAFE} and then to a CALE problem \eqref{eq:CALE}; then we can use globally guaranteed low-rank projected gradient descent (Algorithm \ref{algorithm:LPGD}) to efficiently find the global optimum of the resulting CALE problem, and then we can pull that solution back to the original space. While this approach has proven to be quite powerful in analyzing SDL problems in the present work, it is interesting to note that our double-lifting technique does not immediately apply to %all other related models. For example,
the supervised PCA model proposed by Ritchie et al.  [61] for finding low-dimensional subspace that is also effective for a regression task:
\begin{align}\label{eq:SPCA_regression_W1}
\min_{\W\in \R^{p\times r},\, \W^{T}\W=\I_{r},\,  \Beta\in \R^{1\times r}  } \left[  f\left( \W \begin{bmatrix} \Beta ,\,  \W^{T}  \end{bmatrix} \right):= \lVert \Y_{\textup{label}} -  \Beta^{T}  \W^{T}\X_{\textup{data}} \rVert_{F}^{2} + \xi  \lVert \X_{\textup{data}} - \W \W^{T} \X_{\textup{data}}   \rVert_{F}^{2}\right].
\end{align}
Even though we can realize the objective function in the right hand side of \eqref{eq:SPCA_regression_W1} as a function depending on the product $\W[\Beta,\W^{T}]$, the two matrix factors $\W$ and $[\Beta,\W^{T}]$ are \textit{not} decoupled as before. It is for a future investigation to devise a lifting technique for SPCA and related problems, and obtain a strong global convergence guarantee.

\section*{Acknowledgements}

HL is partially supported by NSF  DMS-2206296 and DMS-2010035.

\small{
\bibliographystyle{amsalpha}
\newcommand{\etalchar}[1]{$^{#1}$}
\providecommand{\bysame}{\leavevmode\hbox to3em{\hrulefill}\thinspace}
\providecommand{\MR}{\relax\ifhmode\unskip\space\fi MR }
% \MRhref is called by the amsart/book/proc definition of \MR.
\providecommand{\MRhref}[2]{%
\href{http://www.ams.org/mathscinet-getitem?mr=#1}{#2}
\providecommand{\href}[2]{#2}
\begin{thebibliography}{RWRY11}

[1]
Woody Austin, Dylan Anderson, and Joydeep Ghosh, \emph{Fully supervised
non-negative matrix factorization for feature extraction}, IGARSS 2018-2018
IEEE International Geoscience and Remote Sensing Symposium, IEEE, 2018,

[2]
Alekh Agarwal, Sahand Negahban, and Martin~J Wainwright, \emph{Fast global
convergence rates of gradient methods for high-dimensional statistical
recovery}, Advances in Neural Information Processing Systems \textbf{23}

[3]
Sanae Amani and Christos Thrampoulidis, \emph{Ucb-based algorithms for
multinomial logistic regression bandits}, Advances in Neural Information
Processing Systems \textbf{34} (2021).

[4]
Herv{\'e} Abdi and Lynne~J Williams, \emph{Principal component analysis}, Wiley
interdisciplinary reviews: computational statistics \textbf{2} (2010), no.~4,

[5]
Christopher~M Bishop et~al., \emph{Neural networks for pattern recognition},
Oxford university press, 1995.

[6]
Shivam Bansa, \emph{fake job postings dataset},
\url{https://www.kaggle.com/datasets/shivamb/real-or-fake-fake-jobposting-prediction}

[7]
\bysame, \emph{Fbi report on fake job postings},
\url{https://www.thesslstore.com/blog/fake-jobs-cybercriminals-prey-on-job-seekers-via-fake-job-postings/}

[8]
Michael~W Berry and Murray Browne, \emph{Email surveillance using non-negative
matrix factorization}, Computational \& Mathematical Organization Theory
\textbf{11} (2005), no.~3, 249--264.

[9]
Michael~W Berry, Murray Browne, Amy~N Langville, V~Paul Pauca, and Robert~J
Plemmons, \emph{Algorithms and applications for approximate nonnegative
matrix factorization}, Computational statistics \& data analysis \textbf{52}
(2007), no.~1, 155--173.

[10]
Amir Beck, \emph{First-order methods in optimization}, SIAM, 2017.

[11]
Dimitri~P Bertsekas, \emph{Nonlinear programming}, Journal of the Operational
Research Society \textbf{48} (1997), no.~3, 334--334.

[12]
\bysame, \emph{Nonlinear programming}, Athena scientific Belmont, 1999.

[13]
Rostyslav Boutchko, Debasis Mitra, Suzanne~L Baker, William~J Jagust, and
Grant~T Gullberg, \emph{Clustering-initiated factor analysis application for
tissue classification in dynamic brain positron emission tomography}, Journal
of Cerebral Blood Flow \& Metabolism \textbf{35} (2015), no.~7, 1104--1111.

[14]
Christopher~M Bishop and Nasser~M Nasrabadi, \emph{Pattern recognition and
machine learning}, vol.~4, Springer, 2006.

[15]
David~M Blei, Andrew~Y Ng, and Michael~I Jordan, \emph{Latent dirichlet
allocation}, Journal of machine Learning research \textbf{3} (2003), no.~Jan,

[16]
Dankmar B{\"o}hning, \emph{Multinomial logistic regression algorithm}, Annals
of the institute of Statistical Mathematics \textbf{44} (1992), no.~1,

[17]
Moody~T Chu, Robert~E Funderlic, and Robert~J Plemmons, \emph{Structured low
rank approximation}, Linear algebra and its applications \textbf{366} (2003),

[18]
Yang Chen, Xiao Wang, Cong Shi, Eng~Keong Lua, Xiaoming Fu, Beixing Deng, and
Xing Li, \emph{Phoenix: A weight-based network coordinate system using matrix
factorization}, IEEE Transactions on Network and Service Management
\textbf{8} (2011), no.~4, 334--347.

[19]
Rick Durrett, \emph{Probability: theory and examples}, fourth ed., Cambridge
Series in Statistical and Probabilistic Mathematics, Cambridge University
Press, Cambridge, 2010.

[20]
Michael Elad and Michal Aharon, \emph{Image denoising via sparse and redundant
representations over learned dictionaries}, IEEE Transactions on Image
processing \textbf{15} (2006), no.~12, 3736--3745.

[21]
Jianqing Fan and Runze Li, \emph{Variable selection via nonconcave penalized
likelihood and its oracle properties}, Journal of the American statistical
Association \textbf{96} (2001), no.~456, 1348--1360.

[22]
David~G Feingold and Richard~S Varga, \emph{Block diagonally dominant matrices
and generalizations of the gerschgorin circle theorem.}, Pacific Journal of
Mathematics \textbf{12} (1962), no.~4, 1241--1250.

[23]
Mehrdad~J Gangeh, Pouria Fewzee, Ali Ghodsi, Mohamed~S Kamel, and Fakhri
Karray, \emph{Multiview supervised dictionary learning in speech emotion
recognition}, IEEE/ACM Transactions on Audio, Speech, and Language Processing
\textbf{22} (2014), no.~6, 1056--1068.

[24]
Mehrdad~J Gangeh, Ahmed~K Farahat, Ali Ghodsi, and Mohamed~S Kamel,
\emph{Supervised dictionary learning and sparse representation-a review},
arXiv preprint arXiv:1502.05928 (2015).

[25]
Gene~H Golub and Christian Reinsch, \emph{Singular value decomposition and
least squares solutions}, Linear algebra, Springer, 1971, pp.~134--151.

[26]
Luigi Grippof and Marco Sciandrone, \emph{Globally convergent block-coordinate
techniques for unconstrained optimization}, Optimization methods and software
\textbf{10} (1999), no.~4, 587--637.

[27]
Luigi Grippo and Marco Sciandrone, \emph{On the convergence of the block
nonlinear gauss--seidel method under convex constraints}, Operations research
letters \textbf{26} (2000), no.~3, 127--136.

[28]
Roger~A Horn and Charles~R Johnson, \emph{Matrix analysis}, Cambridge
university press, 2012.

[29]
Jamie Haddock, Lara Kassab, Sixian Li, Alona Kryshchenko, Rachel Grotheer,
Elena Sizikova, Chuntian Wang, Thomas Merkh, RWMA Madushani, Miju Ahn,
et~al., \emph{Semi-supervised nmf models for topic modeling in learning
tasks}, arXiv preprint arXiv:2010.07956 (2020).

[30]
Nathan Halko, Per-Gunnar Martinsson, and Joel~A Tropp, \emph{Finding structure
with randomness: Probabilistic algorithms for constructing approximate matrix
decompositions}, SIAM review \textbf{53} (2011), no.~2, 217--288.

[31]
Prateek Jain, Raghu Meka, and Inderjit Dhillon, \emph{Guaranteed rank
minimization via singular value projection}, Advances in Neural Information
Processing Systems \textbf{23} (2010).

[32]
Prateek Jain, Praneeth Netrapalli, and Sujay Sanghavi, \emph{Low-rank matrix
completion using alternating minimization}, Proceedings of the forty-fifth
annual ACM symposium on Theory of computing, 2013, pp.~665--674.

[33]
Hamed Jelodar, Yongli Wang, Chi Yuan, Xia Feng, Xiahui Jiang, Yanchao Li, and
Liang Zhao, \emph{Latent dirichlet allocation (lda) and topic modeling:
models, applications, a survey}, Multimedia Tools and Applications
\textbf{78} (2019), no.~11, 15169--15211.

[34]
Tamara~G Kolda and Brett~W Bader, \emph{Tensor decompositions and
applications}, SIAM review \textbf{51} (2009), no.~3, 455--500.

[35]
Daniel~S Kermany, Michael Goldbaum, Wenjia Cai, Carolina~CS Valentim, Huiying
Liang, Sally~L Baxter, Alex McKeown, Ge~Yang, Xiaokang Wu, Fangbing Yan,
et~al., \emph{Identifying medical diagnoses and treatable diseases by
image-based deep learning}, Cell \textbf{172} (2018), no.~5, 1122--1131.

[36]
Jingu Kim, Yunlong He, and Haesun Park, \emph{Algorithms for nonnegative matrix
and tensor factorizations: A unified view based on block coordinate descent
framework}, Journal of Global Optimization \textbf{58} (2014), no.~2,

[37]
Donghyun Kim, Chanyoung Park, Jinoh Oh, Sungyoung Lee, and Hwanjo Yu,
\emph{Convolutional matrix factorization for document context-aware
recommendation}, Proceedings of the 10th ACM conference on recommender
systems, 2016, pp.~233--240.

[38]
Koulik Khamaru and Martin Wainwright, \emph{Convergence guarantees for a class
of non-convex and non-smooth optimization problems}, International Conference
on Machine Learning, PMLR, 2018, pp.~2601--2610.

[39]
Guillaume Lecu{\'e} and Shahar Mendelson, \emph{Sparse recovery under weak
moment assumptions}, Journal of the European Mathematical Society \textbf{19}
(2017), no.~3, 881--904.

[40]
Daniel~D Lee and H~Sebastian Seung, \emph{Learning the parts of objects by
non-negative matrix factorization}, Nature \textbf{401} (1999), no.~6755,

[41]
Daniel Lee and H~Sebastian Seung, \emph{Algorithms for non-negative matrix
factorization}, Advances in neural information processing systems \textbf{13}
(2000), 556--562.

[42]
Daniel~D Lee and H~Sebastian Seung, \emph{Algorithms for non-negative matrix
factorization}, Advances in neural information processing systems, 2001,

[43]
Johannes Leuschner, Maximilian Schmidt, Pascal Fernsel, Delf Lachmund, Tobias
Boskamp, and Peter Maass, \emph{Supervised non-negative matrix factorization
methods for maldi imaging applications}, Bioinformatics \textbf{35} (2019),
no.~11, 1940--1947.

[44]
Hanbaek Lyu, \emph{Convergence and complexity of block coordinate descent with
diminishing radius for nonconvex optimization}, arXiv preprint
arXiv:2012.03503 (2020).

[45]
Julien Mairal, \emph{Optimization with first-order surrogate functions},
International Conference on Machine Learning, 2013, pp.~783--791.

[46]
\bysame, \emph{Stochastic majorization-minimization algorithms for large-scale
optimization}, Advances in Neural Information Processing Systems, 2013,

[47]
Jon Mcauliffe and David Blei, \emph{Supervised topic models}, Advances in
neural information processing systems \textbf{20} (2007).

[48]
Julien Mairal, Francis Bach, and Jean Ponce, \emph{Task-driven dictionary
learning}, IEEE transactions on pattern analysis and machine intelligence
\textbf{34} (2011), no.~4, 791--804.

[49]
Julien Mairal, Francis Bach, Jean Ponce, and Guillermo Sapiro, \emph{Online
learning for matrix factorization and sparse coding}, Journal of Machine
Learning Research \textbf{11} (2010), no.~Jan, 19--60.

[50]
Julien Mairal, Michael Elad, and Guillermo Sapiro, \emph{Sparse representation
for color image restoration}, IEEE Transactions on Image Processing
\textbf{17} (2007), no.~1, 53--69.

[51]
Raghu Meka, Prateek Jain, and Inderjit~S Dhillon, \emph{Guaranteed rank
minimization via singular value projection}, arXiv preprint arXiv:0909.5457

[52]
Julien Mairal, Jean Ponce, Guillermo Sapiro, Andrew Zisserman, and Francis
Bach, \emph{Supervised dictionary learning}, Advances in Neural Information
Processing Systems \textbf{21} (2008), 1033--1040.

[53]
Yu~Nesterov, \emph{Gradient methods for minimizing composite functions},
Mathematical programming \textbf{140} (2013), no.~1, 125--161.

[54]
Sahand Negahban and Martin~J Wainwright, \emph{Estimation of (near) low-rank
matrices with noise and high-dimensional scaling}, The Annals of Statistics
\textbf{39} (2011), no.~2, 1069--1097.

[55]
Gabriel Peyr{\'e}, \emph{Sparse modeling of textures}, Journal of Mathematical
Imaging and Vision \textbf{34} (2009), no.~1, 17--31.

[56]
Dohyung Park, Anastasios Kyrillidis, Srinadh Bhojanapalli, Constantine
Caramanis, and Sujay Sanghavi, \emph{Provable non-convex projected gradient
descent for a class of constrained matrix optimization problems}, stat
\textbf{1050} (2016), 4.

[57]
Dohyung Park, Anastasios Kyrillidis, Constantine Carmanis, and Sujay Sanghavi,
\emph{Non-square matrix sensing without spurious local minima via the
burer-monteiro approach}, Artificial Intelligence and Statistics, PMLR, 2017,

[58]
Dohyung Park, Anastasios Kyrillidis, Constantine Caramanis, and Sujay Sanghavi,
\emph{Finding low-rank solutions via nonconvex matrix factorization,
efficiently and provably}, SIAM Journal on Imaging Sciences \textbf{11}
(2018), no.~4, 2165--2204.

[59]
Michael~JD Powell, \emph{On search directions for minimization algorithms},
Mathematical programming \textbf{4} (1973), no.~1, 193--201.

[60]
Fabian Pedregosa, Ga{\"e}l Varoquaux, Alexandre Gramfort, Vincent Michel,
Bertrand Thirion, Olivier Grisel, Mathieu Blondel, Peter Prettenhofer, Ron
Weiss, Vincent Dubourg, et~al., \emph{Scikit-learn: Machine learning in
python}, the Journal of machine Learning research \textbf{12} (2011),

[61]
Alexander Ritchie, Laura Balzano, Daniel Kessler, Chandra~S Sripada, and
Clayton Scott, \emph{Supervised pca: A multiobjective approach}, arXiv
preprint arXiv:2011.05309 (2020).

[62]
Benjamin Recht, Maryam Fazel, and Pablo~A Parrilo, \emph{Guaranteed
minimum-rank solutions of linear matrix equations via nuclear norm
minimization}, SIAM review \textbf{52} (2010), no.~3, 471--501.

[63]
Bin Ren, Laurent Pueyo, Guangtun~Ben Zhu, John Debes, and Gaspard Duch{\^e}ne,
\emph{Non-negative matrix factorization: robust extraction of extended
structures}, The Astrophysical Journal \textbf{852} (2018), no.~2, 104.

[64]
Pradeep Ravikumar, Martin~J Wainwright, Garvesh Raskutti, and Bin Yu,
\emph{High-dimensional covariance estimation by minimizing ℓ1-penalized
log-determinant divergence}, Electronic Journal of Statistics \textbf{5}
(2011), 935--980.

[65]
Mark Steyvers and Tom Griffiths, \emph{Probabilistic topic models}, Handbook of
latent semantic analysis \textbf{427} (2007), no.~7, 424--440.

[66]
Arkadiusz Sitek, Grant~T Gullberg, and Ronald~H Huesman, \emph{Correction for
ambiguous solutions in factor analysis using a penalized least squares
objective}, IEEE transactions on medical imaging \textbf{21} (2002), no.~3,

[67]
Tao Sun, Yuejiao Sun, and Wotao Yin, \emph{On markov chain gradient descent},
Advances in Neural Information Processing Systems, 2018, pp.~9896--9905.

[68]
Stephen Tu, Ross Boczar, Max Simchowitz, Mahdi Soltanolkotabi, and Ben Recht,
\emph{Low-rank solutions of linear matrix equations via procrustes flow},
International Conference on Machine Learning, PMLR, 2016, pp.~964--973.

[69]
Gui-Xian Tian and Ting-Zhu Huang, \emph{Inequalities for the minimum eigenvalue
of m-matrices}, The Electronic Journal of Linear Algebra \textbf{20} (2010),

[70]
Tian Tong, Cong Ma, and Yuejie Chi, \emph{Accelerating ill-conditioned low-rank
matrix estimation via scaled gradient descent}, Journal of Machine Learning
Research \textbf{22} (2021), no.~150, 1--63.

[71]
Leo Taslaman and Bj{\"o}rn Nilsson, \emph{A framework for regularized
non-negative matrix factorization, with application to the analysis of gene
expression data}, PloS one \textbf{7} (2012), no.~11, e46331.

[72]
Roman Vershynin, \emph{High-dimensional probability: An introduction with
applications in data science}, vol.~47, Cambridge university press, 2018.

[73]
\emph{{WebMD} coronavirus and pneumonia},
\url{https://www.webmd.com/lung/covid-and-pneumonia#1}, Accessed: 2021-03-13.

[74]
Stephen~J Wright, \emph{Coordinate descent algorithms}, Mathematical
Programming \textbf{151} (2015), no.~1, 3--34.

[75]
Michael~E Wall, Andreas Rechtsteiner, and Luis~M Rocha, \emph{Singular value
decomposition and principal component analysis}, A practical approach to
microarray data analysis, Springer, 2003, pp.~91--109.

[76]
Rachel Ward, Xiaoxia Wu, and Leon Bottou, \emph{Adagrad stepsizes: Sharp
convergence over nonconvex landscapes}, International Conference on Machine
Learning, PMLR, 2019, pp.~6677--6686.

[77]
Lingxiao Wang, Xiao Zhang, and Quanquan Gu, \emph{A unified computational and
statistical framework for nonconvex low-rank matrix estimation}, Artificial
Intelligence and Statistics, PMLR, 2017, pp.~981--990.

[78]
Yi~Xu, Zhuoning Yuan, Sen Yang, Rong Jin, and Tianbao Yang, \emph{On the
convergence of (stochastic) gradient descent with extrapolation for
non-convex optimization}, arXiv preprint arXiv:1901.10682 (2019).

[79]
Pavel Yaskov, \emph{Controlling the least eigenvalue of a random gram matrix},
Linear Algebra and its Applications \textbf{504} (2016), 108--123.

[80]
Yael Yankelevsky and Michael Elad, \emph{Structure-aware classification using
supervised dictionary learning}, 2017 IEEE International Conference on
Acoustics, Speech and Signal Processing (ICASSP), IEEE, 2017, pp.~4421--4425.

[81]
Shijie Zhao, Junwei Han, Jinglei Lv, Xi~Jiang, Xintao Hu, Yu~Zhao, Bao Ge, Lei
Guo, and Tianming Liu, \emph{Supervised dictionary learning for inferring
concurrent brain networks}, IEEE transactions on medical imaging \textbf{34}
(2015), no.~10, 2036--2045.

[82]
Qiang Zhang and Baoxin Li, \emph{Discriminative k-svd for dictionary learning
in face recognition}, 2010 IEEE computer society conference on computer
vision and pattern recognition, IEEE, 2010, pp.~2691--2698.

[83]
Qinqing Zheng and John Lafferty, \emph{A convergent gradient descent algorithm
for rank minimization and semidefinite programming from random linear
measurements}, arXiv preprint arXiv:1506.06081 (2015).

[84]
Tuo Zhao, Zhaoran Wang, and Han Liu, \emph{A nonconvex optimization framework
for low rank matrix estimation}, Advances in Neural Information Processing
Systems \textbf{28} (2015), 559.

\end{thebibliography}

\addresseshere

\newpage

\appendix

\section{Proof of main results}

\subsection{Proof of Theorem \ref{thm:CALE_LPGD}}
We first establish Theorem \ref{thm:CALE_LPGD}, which shows exponential convergence  of the low-rank projected gradient descent (Algorithm \ref{algorithm:LPGD}) for the CALE problem \ref{eq:CALE}. The proof is similar to the standard argument that shows exponential convergence  projected gradient descent with fixed step size for constrained strongly convex problems (see, e.g., \cite[Thm. 10.29]{beck2017first}). However, when we have a strongly convex minimization problem with low-rank-constrained matrix parameter, then the constraint set of low-rank matrices is not convex, so one cannot use non-expansiveness of convex projection operator. Indeed, the rank-$r$ projection $\Pi_{r}$ by truncated SVD is not guaranteed to be non-expansive. In order to circumvent this issue, we use the idea of approximating rank-$r$ projection by a suitable linear projection on a carefully chosen linear subspace, an approach used in [77]. Then one can show that rank-$r$ projection is at most $2$-Lipschitz in some sense, so if the contraction constant in standard analysis of projected gradient descent for strongly convex objectives is small enough ($<1/2$), then overall one still retains exponential convergence.

\begin{lemma}(Linear approximation of rank-$r$ projection)\label{lem:rank_r_lin_appx}
Fix $\Y\in \R^{d_{1}\times d_{2}}$, $R\ge r \in \mathbb{N}$, and denote $\X=\Pi_{r}(\Y)$ and $\hat{\X} = \Pi_{\mathcal{A}}(
\Y)$, where $\mathcal{A}\subseteq \R^{d_{1}\times d_{2}}$ is a linear subspace. Let $\X=\U\bSigma \V^{T}$ denote the SVD of $\X$. Suppose there exists $\overline{\U}\in \R^{d_{1}\times R}$ and $\overline{\V} \in \R^{d_{2}\times R}$ such that
\begin{align}
\mathcal{A} = \left\{ \A \in \R^{d_{1}\times d_{2}}\,\big|\, \textup{col}(\A^{T}) \subseteq \textup{col}(\overline{\V}),\, \textup{col}(\A) \subseteq \textup{col}(\overline{\U})  \right\}, \quad \textup{col}(\U)\subseteq \textup{col}(\overline{\U}), \quad \textup{col}(\V)\subseteq \textup{col}(\overline{\V}).
\end{align}
Then $\X=\Pi_{r}(\hat{\X})$.
\end{lemma}

\begin{proof}
Write $\Y-\X=\dot{\U} \dot{\bSigma} \dot{\V}^{T}$  for its SVD. Let $d:=\rank(\Y)$ and let $\sigma_{1}\ge \dots \ge \sigma_{d}>0$ denote the nonzero singular values of $\Y$. Since $\X=\Pi_{r}(\Y) = \U\bSigma\V^{T}$ and $\Y=\U\bSigma\V^{T} + \dot{\U} \dot{\bSigma} \dot{\V}^{T}$, we must have that $\bSigma$ consists of the top $r$ singular values of $\Y$ and the rest of $d-r$ singular values are contained in $\dot{\bSigma}$. Furthermore, $\textup{col}(\U) \perp \textup{col}(\dot{\U})$.

Now, since $\X\in \mathcal{A}$ and $\Pi_{\mathcal{A}}$ is linear, we get
\begin{align}\label{eq:3r_r_approx1_lem}
\hat{\X} = \Pi_{\mathcal{A}}( \X + (\Y -\X) ) = \U \bSigma \V^{T} + \Pi_{\mathcal{A}}(\dot{\U} \dot{\bSigma} \dot{\V}^{T}).
\end{align}
Let $\bZ:= \Pi_{\mathcal{A}}(\dot{\U} \dot{\bSigma} \dot{\V}^{T})$ and write its SVD as $\bZ = \widetilde{\U}\widetilde{\bSigma}\widetilde{\V}^{T}$. Then note that $(\U^{T}\overline{\U} \, \overline{\U}^{T})^{T}=\overline{\U}\,\overline{\U}^{T}\U=\U$ since $\overline{\U}\,\overline{\U}^{T}:\R^{d_{1}}\rightarrow \R^{d_{1}}$ is the orthogonal projection onto $\textup{col}(\overline{\U})\supseteq \textup{col}(\U)$. Hence $\U^{T}\overline{\U}\,\overline{\U}^{T} = \U^{T}$, so we get
\begin{align}
\U^{T} \bZ =  	\left( \U^{T} \overline{\U}\, \overline{\U}^{T}\right) \dot{\U} \dot{\bSigma} \dot{\V}^{T} 	\V^{T} \overline{\V} = \left( \U^{T} \dot{\U} \right) \dot{\bSigma} \dot{\V}^{T} 	\V^{T} \overline{\V} = O.
\end{align}
It follows that $\U^{T}\widetilde{\U}=O$, since $	\U^{T}\widetilde{\U} = \U^{T} \bZ \widetilde{\V} (\widetilde{\bSigma})^{-1} = O$. Therefore, rewriting \eqref{eq:3r_r_approx1_lem} gives the SVD of $\hat{\X}$ as
\begin{align}
\hat{\X} = \begin{bmatrix}
\U & \widetilde{\U}
\end{bmatrix}
\begin{bmatrix}
\bSigma & O\\
O & \widetilde{\bSigma}
\end{bmatrix}
\begin{bmatrix}
\V \\ \widetilde{\V}
\end{bmatrix}
\end{align}
Furthermore, $\lVert  \Pi_{\mathcal{A}}(\dot{\U} \dot{\bSigma} \dot{\V}^{T}) \rVert_{2} \le \lVert\dot{\bSigma} \rVert_{2} =\sigma_{r+1}^{t}$, so $\Sigma$ consists of the top $r$ singular values of $\hat{\X}$. It follows that $\X=\U \bSigma \V^{T}$ is the best rank-$r$ approximation of $\hat{\X}$, as desired.
\end{proof}

\begin{proof}[\textbf{Proof of Theorem} \ref{thm:CALE_LPGD}]
Denote $\bZ^{\star}=[\X^{\star}, \bGamma^{\star}]\in \Param\subseteq \R^{d_{1}\times d_{2}}\times \R^{d_{3}\times d_{4}}$.  Let $\mathcal{A}$ denote a linear subspace of $\R^{d_{1}\times d_{2}}$. Denote
\begin{align}
\hat{\bZ}_{t} =\Pi_{\mathcal{A}\times \R^{d_{3}\times d_{4}}}\left(  \Pi_{\Param} \left( \bZ_{t-1} -  \tau \nabla f(\bZ_{t}) \right) \right).
\end{align}
We will choose $\mathcal{A}$ in such a way that
\begin{align}\label{eq:mathcal_A_subspace_cond}
\X_{t} = \Pi_{r}(\hat{\X}_{t}) \in \argmin_{\X, \rank(\X)\le r} \lVert \hat{\X}_{t} - \X  \rVert_{F}, \quad  \bZ^{\star}\in \mathcal{A} \times \R^{d_{3}\times d_{4}}.
\end{align}
For instance, $\mathcal{A}=\R^{d_{1}\times d_{2}}$ satisfies the above conditions, although this choice does not give optimal control on the variance term (the second term in the right hand side of \eqref{eq:CALE_LPGD_thm}). We will first derive a general bound with such $\mathcal{A}$, and at the end of the proof, we will give a specific construction of such $\mathcal{A}$ to obtain the bound in the assertion.

Denote $\Delta \bZ^{\star}:=\bZ^{\star} - \Pi_{\Param}\left( \bZ^{\star} - \tau \nabla f(\bZ^{\star}))\right)$. Using $\bZ^{\star}\in \mathcal{A}\times \R^{d_{3}\times d_{4}}$ and linearity of the linear projection $\Pi_{\mathcal{A}\times \R^{d_{3}\times d_{4}}}$, write
\begin{align}
\bZ^{\star} &= \Pi_{\mathcal{A}\times \R^{d_{3}\times d_{4}}}(\bZ^{\star}) \\
&=\Pi_{\mathcal{A}\times \R^{d_{3}\times d_{4}}} \left( \Pi_{\Param}(\bZ^{\star} - \tau \nabla f(\bZ^{\star})) \right)  + \Pi_{\mathcal{A}\times \R^{d_{3}\times d_{4}}} \left( \bZ^{\star}-  \Pi_{\Param}(\bZ^{\star} - \tau \nabla f(\bZ^{\star}))  \right) \\
&=\Pi_{\mathcal{A}\times \R^{d_{3}\times d_{4}}} \left( \Pi_{\Param}(\bZ^{\star} - \tau \nabla f(\bZ^{\star})) \right) +  \Pi_{\mathcal{A}\times \R^{d_{3}\times d_{4}}} \left(  \Delta \bZ^{\star} \right).
\end{align}
Namely, the first term above is a one-step update of a projected gradient descent at $\bZ^{\star}$ over $\Param$ with stepsize $\tau$, and the second term above is the error term. If $\bZ^{\star}$ is a stationary point of $f$ over $\Param$, then $-\nabla f(\bZ^{\star})$ lies in the normal cone of $\Param$ at $\bZ^{\star}$, so $\bZ^{\star}$ is invariant under the projected gradient descent and the error term above (the second term in the last expression) is zero. If $\bZ^{\star}$ is only approximately stationary, then the error above is nonzero.

Now,  recall that $\hat{Z}_{t}$ is obtained by using the orthogonal projection $\Pi_{\mathcal{A}}$ onto the convex subset $\mathcal{A}$ instead of the rank-$r$ projection $\Pi_{r}$  to obtain the matrix coordinate of $\hat{\bZ}_{t}$. Notice that $\Pi_{\Param}$ and $\Pi_{\mathcal{A}\times \R^{d_{3}\times d_{4}}}$ are non-expansive being projection onto a convex set, while the rank-$r$ projection $\Pi_{r}$ is not in general. Also using the linearity of the  subspace projection $\Pi_{\mathcal{A}\times \R^{d_{3}\times d_{4}}}$, we get
\begin{align}
&\lVert \hat{\bZ}_{t} - \bZ^{\star} \rVert_{F}  \\
&\qquad = \left\lVert   \Pi_{\mathcal{A}\times \R^{d_{3}\times d_{4}}}\left(  \Pi_{\Param} \left( \bZ_{t-1} -  \tau \nabla f(\bZ_{t-1}) \right) \right) - \Pi_{\mathcal{A}\times \R^{d_{3}\times d_{4}}}\left(  \Pi_{\Param} \left( \bZ^{\star} -  \tau \nabla f(\bZ^{\star}) \right) \right) +\Pi_{\mathcal{A}\times \R^{d_{3}\times d_{4}}}\left( \Delta \bZ^{\star} \right)  \right\rVert_{F}  \\
&\qquad \le  \left\lVert  \bZ_{t-1} -  \tau \nabla f(\bZ_{t-1})  - \bZ^{\star} +  \tau \nabla f(\bZ^{\star})  \right\rVert_{F}  + \lVert \Pi_{\mathcal{A}\times \R^{d_{3}\times d_{4}}}\left( \Delta \bZ^{\star}  \right)  \rVert_{F} \\
&\qquad \le \max( |1-\tau L| ,\, |1-\tau \mu|  ) \, \lVert  \bZ_{t-1} - \bZ^{\star}\rVert_{F} + \lVert \Pi_{\mathcal{A}\times \R^{d_{3}\times d_{4}}}\left( \Delta \bZ^{\star} \right)  \rVert_{F}.
\end{align}

The last inequality follows from the fact that $\bZ_{t}$ and $\bZ^{\star}$ have rank $\le r$ and the restricted strong convexity and smoothness properties (Definition \ref{def:RSC}). Namely, fix $\X,\Y\in \R^{d_{1}\times d_{2}}\times \R^{d_{3}\times d_{4}}$ whose first matrix components have rank $\le r$. Assuming $\nabla f$ is continuous,
\begin{align}
\X -  \tau \nabla f(\X)  - \Y +  \tau \nabla f(\Y) &= (\X-\Y) - \tau ( \nabla f(\X) -  \nabla f(\Y)) \\
&= \int_{0}^{1} \left(  \I-\tau \nabla^{2}(\X + s(\Y-\X)) \right)(\X-\Y)\,ds.
\end{align}
Using the inequality $\lVert \A \B \rVert_{F} \le \lVert \A \rVert_{2} \lVert \B \rVert_{F}$, this gives
\begin{align}
\lVert \X -  \tau \nabla f(\X)  - \Y +  \tau \nabla f(\Y) \rVert_{F} & \le \sup_{\bZ=[\bZ_{1},\bZ_{2}]: \, \rank(\bZ_{1})\le r}   \lVert  \I-\tau \nabla^{2}f( \bZ)\rVert_{2}  \, \lVert \X - \Y \rVert_{F} \\
&\le \rho \,\lVert \X - \Y \rVert_{F},
\end{align}
where $\eta:=\max( |1-\tau L| ,\, |1-\tau \mu| )$. Indeed, for the second inequality above, note that the eigenvalues of $\nabla^{2} f(\bZ)$ are contained in $[\mu,L]$, so the eigenvalues of $\I-\tau \nabla^{2}f(\bZ) $ are between $\min(1-\tau L,\, 1-\tau \mu)$ and $\max(1-\tau L,\, 1-\tau \mu)$. Combining with the previous inequality, it follows that
\begin{align}\label{eq:pf_LPGD_pf1}
\lVert \hat{\bZ}_{t} - \bZ^{\star}  \rVert_{F}  \le  \eta \, \lVert  \bZ_{t-1} - \bZ^{\star}\rVert_{F} + \lVert \Pi_{\mathcal{A}\times \R^{d_{3}\times d_{4}}}\left(\Delta \bZ^{\star}\right)  \rVert_{F}.
\end{align}
Next, recall the notation $\bZ_{t}=[\X_{t}, \bGamma_{t}]$ and $\hat{\bZ}_{t}=[\hat{\X}_{t}, \bGamma_{t}]$.  By construction, we have $\X_{t} =\Pi_{r}( \hat{\X}_{t})$, so $\X_{t}$ is the best rank-$r$ approximation of $\hat{\X}_{t}$ in the sense that  $\X_{t}=\argmin_{\X, \rank(\X)\le r} \lVert \hat{\X}_{t} - \X  \rVert_{F}$. Then observe that
\begin{align}
\lVert \bZ_{t} - \bZ^{\star}  \rVert_{F}  &\le 	\lVert \bZ_{t} - \hat{\bZ}_{t}  \rVert_{F} + 	\lVert \hat{\bZ}_{t} - \bZ^{\star}  \rVert_{F}  \\
&= 	\lVert \X_{t} - \hat{\X}_{t}  \rVert_{F} + 	\lVert \hat{\bZ}_{t} - \bZ^{\star}  \rVert_{F}  \\
&\le 	\lVert\X^{\star} -  \hat{\X}_{t}  \rVert_{F}  + 	\lVert \hat{\bZ}_{t} - \bZ^{\star}  \rVert_{F}   \le 2 	\lVert \hat{\bZ}_{t} - \bZ^{\star}  \rVert_{F},
\end{align}
so by combining with \eqref{eq:pf_LPGD_pf1}, we get
\begin{align}
\lVert \bZ_{t} - \bZ^{\star}  \rVert_{F}  \le  2\eta \, \lVert  \bZ_{t-1} - \bZ^{\star}\rVert_{F}+ \lVert \Pi_{\mathcal{A}\times \R^{d_{3}\times d_{4}}}\left(\Delta \bZ^{\star} \right)  \rVert_{F}.
\end{align}
Note that $0\le \eta<1/2$ if and only if $\tau\in (\frac{1}{2\mu}, \frac{3}{2L})$, and this interval is non-empty if and only if $L/\mu<3$. Hence for such choice of $\tau$, $0<2\eta<1$,  so
by a recursive application of the above inequality, we obtain
\begin{align}\label{eq:linear_conv_ineq_pf1}
\lVert \bZ_{t} - \bZ^{\star}  \rVert_{F}  \le  (2\eta)^{t}\, \lVert  \bZ_{0} - \bZ^{\star}\rVert_{F} +  \frac{1}{1-2\eta} \lVert \Pi_{\mathcal{A}\times \R^{d_{3}\times d_{4}}}\left(\Delta \bZ^{\star} \right)  \rVert_{F}.
\end{align}

Finally, we bound the variance term in the last expression by choosing a suitable linear subspace $\mathcal{A}\subseteq \R^{d_{1}\times d_{2}}$ satisfying \eqref{eq:mathcal_A_subspace_cond}. Note that $\Delta \bZ^{\star}=\tau [\Delta \X^{\star}, \Delta \bGamma^{\star}]$, where the latter is defined in the statement of Theorem \ref{thm:CALE_LPGD}. Recall that $\bZ^{\star}=[\X^{\star}, \bGamma^{\star}]$. Let $\X^{\star} = \U^{\star} \bSigma^{\star} (\V^{\star})^{T}$ denote the SVD of $\X^{\star}$.  For each iteration $t$, denote $\bZ_{t} = [\X_{t}, \bGamma_{t}]$ and let $\X_{t} = \U_{t} \bSigma_{t} \V_{t}^{T}$ denote the SVD of $\X_{t}$. Since $\X_{t}$ and $\X^{\star}$ have rank at most $r$, all of both $\U^{\star}$, $\U_{t}$, $\V^{\star}$, and $\V_{t}$ have at most $r$ columns. Define a matrix $\U_{3r}$ so that its columns form a basis for the subspace spanned by the columns of $[\U^{\star}, \U_{t-1}, \U_{t}]$. Then $\U_{3r}$ has at most $3r$ columns. Similarly, let  $\U_{3r}$ be a matrix so that its columns form a basis for the subspace spanned by the columns of $[\V^{\star}, \V_{t-1}, \V_{t}]$. Then $\V_{3r}$ has at most $3r$ columns. Now, define the subspace
\begin{align}\label{eq:mathcal_A_construction}
\mathcal{A} := \left\{ \Delta\in \R^{d_{1}\times d_{2}}\,|\, \textup{span}(\Delta^{T}) \subseteq \textup{span}(\V_{3r}),\, \textup{span}(\Delta) \subseteq \textup{span}(\U_{3r})  \right\}.
\end{align}
Note that $\mathcal{A}$ is a convex subset of $\R^{d_{1}\times d_{2}}$.  Also note that, by definition, $\X^{\star}, \X_{t}, \X_{t-1}\in \mathcal{A}$. Let $\Pi_{\mathcal{A}}$ denote the projection operator onto $\mathcal{A}$. More precisely, for each $\X\in \R^{d_{1}\times d_{2}}$, we have
\begin{align}
\Pi_{\mathcal{A}}(\X) = \U_{3r}\U_{3r}^{T} \X \V_{3r}\V_{3r}^{T}.
\end{align}
Then by Lemma \ref{lem:rank_r_lin_appx}, we have $\X_{t} =\Pi_{r}( \hat{\X}_{t})$. Hence $\mathcal{A}$  in \eqref{eq:mathcal_A_construction} satisfies \eqref{eq:mathcal_A_subspace_cond}. Therefore, \eqref{eq:linear_conv_ineq_pf1} holds for the $\mathcal{A}$ chosen as in \eqref{eq:mathcal_A_construction}.

Now, note that  $\Pi_{\mathcal{A}\times \R^{d_{3}\times d_{4}}}(\Delta \X^{\star}, \Delta \bGamma^{\star}) = [\Pi_{\mathcal{A}}(\Delta \X^{\star}) , \Delta \bGamma^{\star} ]$ and  $\rank(\mathcal{A})\le 3r$. Thus by triangle inequality,
\begin{align}\label{eq:LPGD_last_gap_bound}
\lVert \Pi_{\mathcal{A}\times \R^{d_{3}\times d_{4}}}\left( \Delta \X^{\star}, \Delta \bGamma^{\star}  \right)  \rVert_{F} &\le \lVert \Pi_{\mathcal{A}}( \Delta \X^{\star})  \rVert_{F} + \lVert  \Delta \bGamma^{\star} \rVert_{F} \le \sqrt{ 3r} \lVert \Delta \X^{\star} \rVert_{2} +  \lVert \Delta \bGamma^{\star} \rVert_{F}.
\end{align}
This completes the proof of \textbf{(i)}.

Next, we show \textbf{(ii)}. Suppose $\mathbf{Z}^{\star}$ is a stationary point of $f$ over $\Param$. Then $\Delta \bZ^{\star}=O$ so the first part of the assertion follows from \textbf{(i)}. For the second part, suppose that $\nabla f$ is $L'$-Lipschiz over $\Param$ for some $L'>0$. Then by Cauchy-Schwarz inequality,
\begin{align}
\left| f(\bZ_{n}) - f(\bZ^{\star}) \right| & = \left|  \int_{0}^{1} \left\langle \nabla f \left(\bZ_{n} + s(\bZ^{\star} - \bZ_{n}) \right),\, \bZ_{n}-\bZ^{\star} \right\rangle \,ds \right| \\
&\le   \int_{0}^{1}\left\lVert  \nabla f \left(\bZ_{n} + s(\bZ^{\star} - \bZ_{n}) \right) \right\rVert \lVert \bZ_{n}-\bZ^{\star} \rVert \,ds \\
&\le   \int_{0}^{1} \left( \left\lVert  \nabla f (\bZ^{\star}) \rVert + s L\lVert \bZ_{n}-\bZ^{\star} \right\rVert  \right) \lVert \bZ_{n}-\bZ^{\star} \rVert \,ds \\
&\le   \left( \left\lVert  \nabla f (\bZ^{\star}) \rVert +  L\lVert \bZ_{n}-\bZ^{\star} \right\rVert  \right) \lVert \bZ_{n}-\bZ^{\star} \rVert.
\end{align}
Then \eqref{eq:PSGD_linear_conv2} follows by combining the above inequality with \textbf{(i)}.
\end{proof}

\begin{remark}\label{rmk:pf_thm_LPGD}
Note that in \eqref{eq:LPGD_last_gap_bound}, we could have used the following crude bound
\begin{align}\label{eq:LPGD_last_gap_bound2}
\left\lVert \Pi_{\mathcal{A}\times \R^{d_{3}\times d_{4}}}\left( \Delta \X^{\star}, \Delta \bGamma^{\star}  \right)  \right\rVert_{F} \le 	\left\lVert \left[  \Delta \X^{\star}, \Delta \bGamma^{\star}  \right]  \right\rVert_{F} &\le
\lVert \Delta \X^{\star} \rVert_{F} +  \lVert \Delta \bGamma^{\star} \rVert_{F} \\
&\le \sqrt{\rank(\Delta \X^{\star})}	\lVert \Delta \X^{\star} \rVert_{2} +  \lVert \Delta \bGamma^{\star} \rVert_{F},
\end{align}
which is also the bound we would have obtained if we choosed the trivial linear subspace $\mathcal{A}=\R^{d_{1}\times d_{2}}$ in the proof of Theorem \ref{thm:CALE_LPGD} above. While we know $\rank(\X^{\star})\le r $, we do not have an a priori bound on $ \rank(\Delta \X^{\star})$, which could be much larger then $\sqrt{3r}$. A smarter choice of the subspace $\mathcal{A}$ as we used in the proof of Theorem \ref{thm:CALE_LPGD} ensures that we only need the factor $\sqrt{3r}$ in place of the unknown factor $\sqrt{\rank(\Delta \X^{\star})}$ as in \eqref{eq:LPGD_last_gap_bound}.
\end{remark}

\subsection{Proof of Theorems \ref{thm:SDL_LPGD} and \ref{thm:SDL_LPGD_feat}}

Next, prove Theorems \ref{thm:SDL_LPGD} and \ref{thm:SDL_LPGD_feat}, which amounts to verify that the hypothesis of Theorem \ref{thm:CALE_LPGD} holds for the SDL problems in \eqref{eq:SDL_filt_1} and \eqref{eq:SDL_feat_1}.

We begin with some preliminary computation. Let $\a_{s}$ denote the activation corresponding to the $s$th sample (see \eqref{eq:ASDL_1}). More precisely,  $\a_{s}=\A^{T}\x_{s}+\bGamma^{T}\x'_{s}$ for the filter-based model with $\A\in \R^{p\times \kappa}$, and $\a_{s}=\A[:,s]+\bGamma^{T}\x'_{s}$ with $\A\in \R^{\kappa\times n}$. In both cases, $\B\in \R^{p\times n}$ and $\bGamma\in \R^{q\times \kappa}$. Then the objective function $f$ in \eqref{eq:SDL_filt_CALE}   can be written as
\begin{align}\label{eq:f_SDL_pf0}
f(\A, \B, \bGamma)&:=\left( -\sum_{s=1}^{n} \sum_{j=0}^{\kappa} \mathbf{1}(y_{i}=j)  \log g_{j}(  \a_{s}  ) \right) +   \xi \lVert  \X_{\textup{data}}  -\B\rVert_{F}^{2} + \nu \left( \lVert \A \rVert_{F}^{2}+\lVert \bGamma \rVert_{F}^{2} \right) \\
&= \sum_{s=1}^{n} \left(  \log \left( 1+\sum_{c=1}^{\kappa} h( \a_{s}[c]  ) \right) -  \sum_{j=1}^{\kappa} \mathbf{1}(y_{i}=j) \log  h( \a_{s}[j]  )  \right) +   \xi \lVert  \X_{\textup{data}}  -\B\rVert_{F}^{2} + \nu \left( \lVert \A \rVert_{F}^{2} + \lVert \bGamma \rVert_{F}^{2}  \right),
\end{align}
where $\a_{s}[i]\in \R$ denotes the $i$th component of $\a_{s}\in \R^{\kappa}$. In the proofs we provided below, we compute the Hessian of $f$ above explicitly for the filter- and the feature-based cases and use Theorem \ref{thm:CALE_LPGD} to derive the result. Recall the functions $\dot{\h}$ and $\ddot{\H}$ introduced in \ref{assumption:A4}. For each label $y\in \{0,\dots,\kappa\}$ and activation $\a\in \R^{\kappa}$, the negative log likelihood of observing label $y$ from the probability distribution $\g(\a)$ defined in \eqref{eq:ell_log_likelihood} can be written as
\begin{align}
\ell_{0}(y,\a) := \log\left( \sum_{c=1}^{\kappa} h(\a[c]) \right)  - \sum_{c=1}^{\kappa}\mathbf{1}(y=c) \log h(\a[c]).
\end{align}
Then we have the following relations
\begin{align}
\nabla_{\a} \ell_{0}(y,\a) = \dot{\h}(y,\a), \qquad \nabla_{\a}\nabla_{\a^{T}} \ell_{0}(y,\a) = \ddot{\H}(y,\a).
\end{align}

\begin{proof}[\textbf{Proof of Theorem} \ref{thm:SDL_LPGD}]
Let $f=f_{\textup{SDL-filt}}$ denote the loss function for the filter-based SDL model in \eqref{eq:SDL_filt_CALE}. Fix $\bZ_{1},\bZ_{2}\in \Param\subseteq \R^{d_{1}\times d_{2}}\times \R^{d_{3}\times d_{4}} $. By \ref{assumption:A1} $\Param$ is convex, so $t\bZ_{1} + (1-t) \bZ_{2}\in \Param$ for all $t\in [0,1]$. Then by the mean value theorem, there exists $t^{*}\in [0,1]$ such that for $\bZ^{*}=t^{*}\bZ_{1} + (1-t^{*})\bZ_{2}$,
\begin{align}\label{eq:RSC_thm_pf}
f(\bZ_{2}) - f(\bZ_{1}) - \langle \nabla f(\bZ_{1}), \, \bZ_{2}-\bZ_{1} \rangle =  \left( \vect(\bZ_{2}) - \vect(\bZ_{1})  \right)^{T} \nabla_{\vect(\bZ)}\nabla_{\vect(\bZ)^{T}} f( \bZ^{*} )  \left( \vect(\bZ_{2}) - \vect(\bZ_{1})  \right).
\end{align}
Hence, according to Theorem \ref{thm:CALE_LPGD}, it suffices to verify that
for some $\mu,L>0$ such that $L/\mu<3$,
\begin{align}
\frac{\mu}{2} \I \preceq \nabla_{\vect(\bZ)}\nabla_{\vect(\bZ)^{T}} f( \bZ^{*} ) \preceq \frac{L}{2} \I
\end{align}
for all $\bZ^{*}=[\X,\bGamma]$ with $\rank(\X^{*})\le r$.

%According to Theorem \ref{thm:CALE_LPGD}. it boils down to computing the gradient and the Hessian of the SDL loss function $f_{\textup{SDL-filt}}$ in \eqref{eq:SDL_filt_CALE}.

To this end, let $\a_{s}$ denote the activation corresponding to the $s$th sample (see \eqref{eq:ASDL_1}). More precisely,  $\a_{s}=\A^{T}\x_{s}+\bGamma^{T}\x'_{s}$ for the filter-based model we consider here. We discussed that the objective function $f$ in \eqref{eq:SDL_filt_CALE} can be written as \eqref{eq:f_SDL_pf0}. Denote
\begin{align}\label{eq:thm_SDL_filt_as_def}
\a_{s}= \A^{T}\x_{s}+\bGamma^{T}\x_{s}' =: \vast[ \left\langle \underbrace{\begin{bmatrix}
\A[:,j] \\ \bGamma[:,j] \end{bmatrix}}_{=:\u_{j}} ,\, \underbrace{\begin{bmatrix}
\x_{s} \\ \x_{s}' \end{bmatrix}}_{=:\bphi_{s}}  \right\rangle; \,\, j=1,\dots,\kappa \vast]^{T}\in \R^{\kappa},
\end{align}
where we have introduced the notations $\u_{j}\in \R^{(p+q)\times 1}$ for $j=1,\dots,\kappa$  and $\bphi_{s}\in \R^{(p+q) \times 1 }$ for $s=1,\dots, n$. Denote $\U:=[\u_{1},\dots,\u_{\kappa}]\in \R^{(p+q)\times \kappa}$, which is a matrix parameter that combines $\A$ and $\bGamma$. Also denote $\bPhi=(\bphi_{1},\dots,\bphi_{n})\in \R^{(p+q)\times n}$ that combined feature matrix of $n$ observations. Then we can compute the gradient and the Hessian of $f$ above as follows:
\begin{align}\label{eq:SDL_filt_gradients}
&\nabla_{\vect(\U)}  f(\U,\B) = \left(  \sum_{s=1}^{n} \dot{\h}(y_{s},\U^{T}\bphi_{s}) \otimes \bphi_{s} \right) + 2\nu \vect(\U), \quad \nabla_{\B}  f(\U,\B) = 2\xi (\B-\X_{\textup{data}}) \\
&\nabla_{\vect(\U)}\nabla_{\vect(\U)^{T}}  f(\U,\B) =  \left( \sum_{s=1}^{n} \ddot{\H}(y_{s},\U^{T}\bphi_{s}) \otimes \bphi_{s}\bphi_{s}^{T}\right) + 2\nu \I_{(p+q)\kappa}, \\
&\nabla_{\vect(\B)}\nabla_{\vect(\B)^{T}}  f(\U,\B) = 2\xi \I_{pn}, \qquad \nabla_{\vect(\B)}\nabla_{\vect(\U)^{T}}  f(\U,\B) = O,
\end{align}
where $\otimes$ above denotes the Kronecker product and the functions $\dot{\h}$ and $\ddot{\H}$ are defined in  \eqref{eq:Hddot_def}.

Recall that the eigenvalues of $\A\otimes \B$, where $\A$ and $\B$ are two square matrices, are given by $\lambda_{i}\mu_{j}$, where $\lambda_{i}$ and $\mu_{j}$ run over all eigenvalues of $\A$ and $\B$, respectively. Hence denoting $\H_{\U}:=\sum_{s=1}^{N} \ddot{\H}(y_{s},\U^{T}\bphi_{s},) \otimes \bphi_{s}\bphi_{s}^{T}$ and using \ref{assumption:A2}-\ref{assumption:A3}, we can deduce
\begin{align}\label{eq:MNL_evals_bounds}
\lambda_{\min}(\H_{\U}) &\ge n \lambda_{\min}\left( n^{-1} \bPhi \bPhi^{T} \right) 	\min_{1\le s \le N,\, \U} \lambda_{\min}\left( \ddot{\H}(y_{s},\bphi_{s}, \U) \right) \ge n  \delta^{-}\alpha^{-} \ge n \mu^{*}>0, \\
\lambda_{\max}(\H_{\U}) &\le n \lambda_{\max}\left( n^{-1}\bPhi \bPhi^{T} \right) 	\max_{1\le s \le N,\, \U} \lambda_{\min}\left( \ddot{\H}(y_{s},\bphi_{s}, \U) \right) \le n \delta^{+}\alpha^{+}\le n L^{*}.
\end{align}
This holds for all $\A,\B,\bGamma$ such that $\rank([\A,\B])\le r$ and under the convex constraint in \eqref{assumption:A1} (also recall that $\U$ is the vertical stack of $\A$ and $\bGamma$). Hence we conclude that the objective function $f_{\textup{SDL-filt}}$ in \eqref{eq:SDL_filt_CALE} verifies RSC and RSM properties (Def. \ref{def:RSC}) with parameters $\mu=\min(2\xi,2\nu+n \mu^{*})$ and $L=\max(2\xi, 2\nu + n L^{*})$. It is straightforward to verify that $L/\mu<3$ if and only if \eqref{eq:thm_SDL_LGPD_cond} holds. This verifies \eqref{eq:RSC_thm_pf} for the chosen parameters $\mu$ and $L$. 	Then the rest follows from Theorem \ref{thm:CALE_LPGD}.
\end{proof}

Next, we prove Theorem \ref{thm:SDL_LPGD_feat}, the exponential convergence of Algorithm \ref{alg:SDL_feat_LPGD} for the feature-based SDL in \eqref{eq:SDL_feat_1}.

\begin{proof}[\textbf{Proof of Theorem} \ref{thm:SDL_LPGD_feat}]
We will use the same setup as in the Proof of Theorem \ref{thm:SDL_LPGD}. The main part of the argument is the computation of the Hessian of loss function $f:=f_{\textup{SDL-feat}}$ in \eqref{eq:SDL_feat_CALE}, which is straightforward but a bit more involved than the corresponding computation for the filter-based case in the proof of Theorem \ref{thm:SDL_LPGD}. To this end, let $\a_{s}$ denote the activation corresponding to the $s$th sample (see \eqref{eq:ASDL_1}). More precisely,  $\a_{s}=\A+\bGamma^{T}\x'_{s}$ for the feature-based model we consider here. Recall the objective function $f$ in  \eqref{eq:SDL_feat_CALE} re-written in \eqref{eq:f_SDL_pf0}. We will compute the gradient and the Hessian of $f$ below.

Recall that for the feature-based model we consider here, we have  $\a_{s}=\A[:,s]+\bGamma^{T}\x'_{s}$, where in this case $\A\in \R^{\kappa\times n}$ (see \eqref{eq:SDL_feat_CALE}). Denote
\begin{align}
\a_{s}=\I_{\kappa} \A[:,s]+\bGamma^{T}\x_{s}' =: \vast[ \left\langle \underbrace{\begin{bmatrix}
\I_{\kappa}[:,j] \\ \bGamma[:,j] \end{bmatrix}}_{=:\v_{j}} ,\, \underbrace{\begin{bmatrix}
\A[:,s] \\ \x_{s}' \end{bmatrix}}_{=:\bpsi_{s}}  \right\rangle; \,\, j=1,\dots,\kappa \vast]^{T}\in \R^{\kappa}.
\end{align}
Note that for the feature-based model here, $\A[:,s]$ is concatenated with the auxiliary covariate $\x'_{s}$, whereas we concatenated $\A[:,j]$ with $\bGamma[:,j]$ for the filter-based case (see \eqref{eq:thm_SDL_filt_as_def})\footnote{This is because for the feature-based model, the column $\A[:,s]\in \R^{\kappa}$ for $s=1,\dots,n$ represent a feature of the $s$th sample, whereas for the filter-based model, $\A[:,j]$ for $j=1,\dots,\kappa$ represents the $j$th filter that is applied to the feature $\x_{s}$ of the $s$th sample.}.

A straightforward computation shows the following gradient formulas:
\begin{align}\label{eq:SDL_feat_gradients}
&\nabla_{\vect(\bGamma)}  f(\A,\B,\bGamma) = \left(  \sum_{s=1}^{n} \dot{\h}(y_{s},\a_{s}) \otimes \x'_{s} \right) + 2\nu \vect(\bGamma), \\
&\nabla_{\vect(\A)}  f(\V,\B) =  \begin{bmatrix}
\dot{\h}(y_{1},\a_{1}) \\
\vdots  \\
\dot{\h}(y_{n},\a_{n})
\end{bmatrix}
+ 2\nu \vect(\A), \quad \nabla_{\B}  f(\A,\B,\bGamma) = 2\xi(\B-\X_{\textup{data}}) \\
&\nabla_{\vect(\bGamma)}\nabla_{\vect(\bGamma)^{T}}  f(\A,\B,\bGamma) =  \left(  \sum_{s=1}^{n} \ddot{\H}(y_{s},\a_{s}) \otimes \x_{s}'(\x_{s}')^{T}\right) + 2\nu \I_{q\kappa}, \\
&\nabla_{\vect(\A)}\nabla_{\vect(\A)^{T}}  f(\A,\B,\bGamma) = \diag\left( \ddot{\H}(y_{1},\a_{1}),\dots, \ddot{\H}(y_{n},\a_{n}) \right)  + 2\nu \I_{\kappa n} \\
&\nabla_{\vect(\bGamma)}\nabla_{\vect(\A)^{T}}  f(\A,\B,\bGamma) =  \left[ \ddot{\H}(y_{1},\a_{1})\otimes \x'_{1}, \dots, \ddot{\H}(y_{1},\a_{n})\otimes \x'_{n} \right] \in \R^{\kappa q \times \kappa n} \\
&\nabla_{\vect(\B)}\nabla_{\vect(\B)^{T}}  f(\A,\B,\bGamma) = 2\xi \I_{pn}, \qquad \nabla_{\vect(\B)}\nabla_{\vect(\V)^{T}}  f(\A,\B,\bGamma) = O.
\end{align}
From this we will compute the eigenvalues of the Hessian $\H_{\textup{feat}}$ of the loss function $f$. In order to illustrate our computation in a simple setting, we first assume $\kappa=1=q$, which corresponds to binary classification $\kappa=1$ with one-dimensional auxiliary covariates $q=1$. In this case, we have
\begin{align}
\H_{\textup{feat}}&:=\nabla_{\vect(\A,\bGamma,\B)} \nabla_{\vect(\A,\bGamma,\B)^{T}} f(\A,\B,\bGamma) \\
\begin{bmatrix}
\ddot{h}(y_{1},\a_{1})+2\nu  & 0  & \dots & 0 &  \ddot{h}(y_{1},\a_{1}) x_{1}'  & O \\
0 & \ddot{h}(y_{2},\a_{2})+2\nu   &\dots & 0 &  \ddot{h}(y_{2},\a_{2}) x_{2}'  & O \\
\vdots & \vdots & \ddots & \vdots & \vdots & \vdots \\
0 & \dots & 0 &  \ddot{h}(y_{n},\a_{n})+2\nu  &  \ddot{h}(y_{n},\a_{n}) x_{n}' & O \\
\ddot{h}(y_{1},\a_{1}) x_{1}' & \ddot{h}(y_{2},\a_{2}) x_{2}' & \dots & \ddot{h}(y_{n},\a_{n}) x_{n}' & \left( \frac{1}{n}\sum_{s=1}^{n}\ddot{h}(y_{s},\a_{s}) (x_{s}')^{2} \right) + 2\nu   & O \\
O&O&\dots &O &O & 2\xi \I_{pn}
\end{bmatrix},
\end{align}
where we denoted $\ddot{h}=\ddot{h}_{11}\in \R$ and $x'_{s}=\x_{s}'\in \R$ for $s=1,\dots,n$. In order to compute the eigenvalues of the above matrix, we will use the following formula for determinant of $3\times 3$ block matrix: ($O$ representing matrices of zero entries with appropriate sizes)
\begin{align}
\det\left( \begin{bmatrix}
A & B & O \\
B^{T} & C & O \\
O & O & D
\end{bmatrix} \right)
= \det\left( C - B^{T}A^{-1}B \right) \det(A) \det(D).
\end{align}
This yields the following simple formula for the characteristic polynomial of $\H_{\textup{feat}}$:
\begin{align}
\det( \H_{\textup{feat}} - \lambda \I) &= \left(  \sum_{s=1}^{n} \ddot{h}(y_{s},\a_{s}) (x_{s}')^{2}  -  \sum_{s=1}^{n} \frac{(\ddot{h}(y_{s},\a_{s}))^{2} (x_{s}')^{2}}{\ddot{h}(y_{s},\a_{s})+2\nu } + 2\nu  - \lambda \right) (2\xi - \lambda) \prod_{s=1}^{n} \left( \ddot{h}(\y_{s},\a_{s}) + 2\nu- \lambda\right)\\
&=\left( \sum_{s=1}^{n}  \frac{2\nu  \ddot{h}(y_{s},\a_{s}) (x_{s}')^{2}}{ \ddot{h}(y_{s},\a_{s})+2\nu } + 2\nu - \lambda \right) (2\xi - \lambda)^{pn} \prod_{s=1}^{n} \left( \ddot{h}(\y_{s},\a_{s}) + 2\nu- \lambda\right).
\end{align}
By \ref{assumption:A4}, we know that $\ddot{h}(y_{s},\a_{s})>0$ for all $s=1,\dots,n$, so the first term in the parenthesis in the above display is lower bounded by $2\nu-\lambda$. It follows that
\begin{align}
\lambda_{\min}(\H_{\textup{feat}}) & \ge \min(2\xi, \alpha^{-}+2\nu), \\
\lambda_{\max}(\H_{\textup{feat}}) &\le \max\left( 2\nu + \alpha^{+}\sum_{s=1}^{n} (x_{s}')^{2} ,\, 2\xi,\, \alpha^{+} + 2\nu\right).
\end{align}

Now we generalize the above computation for general $\kappa,q\ge 1$ case. First note the general form of the Hessian as below:
\begin{align}
&\H_{\textup{feat}}:=\nabla_{\vect(\A,\bGamma,\B)} \nabla_{\vect(\A,\bGamma,\B)^{T}} f(\A,\B,\bGamma) \\
&\quad =
\begin{bmatrix}
\ddot{\H}(y_{1},\a_{1})+2\nu  \I_{\kappa} & 0  & \dots & 0 & (\ddot{\H}(y_{1},\a_{1})\otimes \x_{1}')^{T} & O \\
0 & \ddot{\H}(y_{2},\a_{2})+2\nu  \I_{\kappa}  &\dots & 0 & (\ddot{\H}(y_{2},\a_{2})\otimes \x_{2}')^{T} & O \\
\vdots & \vdots & \ddots & \vdots & \vdots & \vdots \\
0 & \dots & 0 &  \ddot{\H}(y_{n},\a_{n})+2\nu  \I_{\kappa} & (\ddot{\H}(y_{n},\a_{n})\otimes \x_{n}')^{T} & O \\
\ddot{\H}(y_{1},\a_{1})\otimes \x_{1}' & \ddot{\H}(y_{2},\a_{2})\otimes \x_{2}' & \dots & \ddot{\H}(y_{n},\a_{n})\otimes \x_{n}' & \begin{matrix}  \sum_{s=1}^{n}\ddot{\H}(y_{s},\a_{s}) \otimes \x_{s}'(\x_{s})^{T}  \\  + 2\nu  \I_{q\kappa}   \end{matrix} & O \\
O&O&\dots &O &O & 2\xi  \I_{pn}
\end{bmatrix}.
\end{align}
Note that for any square symmetric matrix $B$ and a column vector $\x$ of matching size,
\begin{align}
B\otimes \x \x^{T} - (B\otimes \x)^{T}(B+\lambda \I)^{-1} (B\otimes \x) &= \left( B-B(B+\lambda \I)^{-1} B \right) \otimes (\x\x^{T})  \\
&= (B+\lambda I)^{-1} B \otimes \x\x^{T} \\
& \preceq \I  \otimes \x\x^{T},
\end{align}
where the last diagonal dominace is due to the Woodbury identity for matrix inverse (e.g., see [28]). Hence by a similar computation as before, we obtain
\begin{align}
\det( n\H_{\textup{feat}} - \lambda \I)
%&\quad = \left( \sum_{s=1}^{n} \ddot{\H}(y_{s},\a_{s})\otimes \x_{s}'(\x_{s})^{T}  - \sum_{s=1}^{n}\ddot{\H}(y_{s},\a_{s}) \left( \ddot{\H}(y_{s},\a_{s}) + 2\nu \I\right)^{-1} \ddot{\H}(y_{s},\a_{s}) \otimes \x_{s}'(\x_{s}')^{T}  + 2\nu \I_{q\kappa} - \lambda \right) (2\xi - \lambda) \prod_{s=1}^{n} \left( \ddot{h}(\y_{s},\a_{s}) + 2\nu- \lambda\right)\\
& = \det\left( \sum_{s=1}^{n} 2\nu \left( \ddot{\H}(y_{s},\a_{s}) + 2\nu  \I_{\kappa} \right)^{-1} \ddot{\H}(y_{s},\a_{s}) \otimes \x_{s}'(\x_{s}')^{T}  + (2\nu -\lambda) \I_{q\kappa}  \right) (2\xi n - \lambda)^{pn} \\
&\hspace{5cm} \times \prod_{s=1}^{n} \det\left( \ddot{\H}(\y_{s},\a_{s}) + (2\nu - \lambda)\I_{\kappa} \right).
\end{align}
It follows that
\begin{align}
\lambda_{\min}(\H_{\textup{feat}}) & \ge \min(2\xi,  \alpha^{-}+2\nu), \\
\lambda_{\max}(\H_{\textup{feat}}) &\le \max\left(2\nu +\alpha^{+} n \lambda_{\max}\left( n^{-1}\X_{\textup{aux}} \X_{\textup{aux}}^{T}  \right),\, 2\xi,\,  \alpha^{+} + 2\nu\right).
\end{align}
Then the rest follows from Theorem \ref{thm:SDL_LPGD_feat}.
\end{proof}

\section{Proof of Theorem \ref{thm:SDL_BCD}}

In this section, we prove Theorem \ref{thm:SDL_BCD} only for the case of filter-based SDL in \eqref{eq:ASDL_1}. An almost identical argument will show the assertion for the feature-based case.

Recall the filter-based SDL loss function $f(\A,\B,\bGamma)$ in \eqref{eq:f_SDL_pf0} in terms of the combined variables $[\A,\B,\bGamma]$, where $\A=\W\Beta$ and $\B=\W\H$. For convenience, recall that $\A\in \R^{\kappa\times p}$, $\W\in \R^{p\times r}$, $\Beta\in \R^{r\times \kappa}$, $\H\in \R^{r\times n}$, and $\bGamma\in \R^{q\times n}$.  In the following computations, we will use \textit{commutation matrix} $\C^{(a\times b)}$, which is a special instance of $ab\times ab$ permutation matrix. Namely, for each integers $a,b\ge 1$, there exists a unique matrix $\C^{(a\times b)}\in \{0,1\}^{ab\times ab}$ such that for all $A\in \R^{a\times b}$, we have $\C^{(a,b)}\vect(A) = \vect(A^{T})$. Note that $(\C^{(a,b)})^{T}=\C^{b\times a}$. Furthermore, $(\C^{(a,b)})^{T}\C^{(a,b)}=\I_{ab}$ since $(\C^{(a,b)})^{T}\C^{(a,b)}\vect(A)=\C^{(b,a)}\vect(A^{T})=\vect(A)$. Hence $\C^{(a,b)}$ is positive semi-definite. Throughout this section, we denote $\bZ=[\W,\H,\Beta,\bGamma]$ for the combined SDL parameters.

\begin{lemma}[Derivatives of the filter-based SDL objective in separate variables]
\label{lem:SDL_filt_BCD_derivatives}
Let $L(\bZ)$ denote the objective of the filter-based SDL in \eqref{eq:ASDL_1}. Suppose \ref{assumption:A4} holds. Recall $\dot{\h}$ and $\ddot{\H}$ defined in \eqref{eq:Hddot_def}. Let $\a_{s}:=\Beta^{T}\W^{T}\x_{s}+\bGamma^{T}\x'_{s}$ for $s=1,\dots,n$ and $\K:=[\dot{\h}(y_{1},\a_{1}),\dots, \dot{\h}(y_{1},\a_{1})]\in \R^{\kappa\times n}$. Then we have
\begin{align}\label{eq:SDL_filt_BCD_derivatives1}
\nabla_{\W} \,  L(\bZ) &= \X_{\textup{data}} \K^{T} \Beta^{T} +  2\xi(\W\H-\X_{\textup{data}})\H^{T},\qquad
\nabla_{\Beta} \, L(\bZ)  = \W^{T} \X_{\textup{data}}  \K^{T} \\
\nabla_{\bGamma} \, L(\bZ) &=  \X_{\textup{aux}}\K^{T}, \qquad  \nabla_{\H} \, L(\bZ) =  2\xi \W^{T}(\W\H-\X_{\textup{data}}).
\end{align}
Furthermore, for diagonal terms in the Hessian, we have
\begin{align}
\nabla_{\vect(\W)}\nabla_{\vect(\W)^{T}} \, L(\bZ)
&= (\Beta \otimes \X_{\textup{data}}) \diag\left( \ddot{\H}(y_{1},\a_{1}),\dots,\ddot{\H}(y_{n},\a_{n})\right) \C^{(n,\kappa)} (\Beta \otimes \X_{\textup{data}})^{T} +2\xi (\H\H^{T} \otimes \I_{p} ), \\
\nabla_{\vect(\H)}\nabla_{\vect(\H)^{T}} \, L(\bZ) & = 2\xi (\I_{n}\otimes \W^{T}\W), \\
\nabla_{\vect(\Beta)}\nabla_{\vect(\Beta)^{T}} \, L(\bZ)
&=(\I_{\kappa}\otimes \W^{T}\X_{\textup{data}}) \diag\left( \ddot{\H}(y_{1},\a_{1}),\dots,\ddot{\H}(y_{n},\a_{n})\right)   (\I_{\kappa}\otimes \W^{T}\X_{\textup{data}})^{T}, \\
\nabla_{\vect(\bGamma)}\nabla_{\vect(\bGamma)^{T}} \, L(\bZ)
&= (\I_{r}\otimes \X_{\textup{aux}}) \diag\left( \ddot{\H}(y_{1},\a_{1}),\dots,\ddot{\H}(y_{n},\a_{n})\right)  (\I_{r}\otimes \X_{\textup{aux}})^{T}.
\end{align}
Lastly, for the off-diagonal terms in the Hessian, we have
\begin{align}
\nabla_{\vect(\Beta)}\nabla_{\vect(\W)^{T}} \, L(\bZ)
&= \C^{(\kappa,  n)} (\I_{\kappa}\otimes \W^{T}\X_{\textup{data}})\diag\left( \ddot{\H}(y_{1},\a_{1}),\dots,\ddot{\H}(y_{n},\a_{n})\right)    (\Beta \otimes \X_{\textup{data}})^{T}  +  \C^{(\kappa, r)} (\I_{r}\otimes \X_{\textup{data}} \K^{T} )^{T}, \\
\nabla_{\vect(\bGamma)}\nabla_{\vect(\W)^{T}} \, L(\bZ) &= O,\\
\nabla_{\vect(\H)}\nabla_{\vect(\W)^{T}} \, L(\bZ) &= 2\xi \left[ (\H^{T}\otimes \W^{T}) + (\I_{r}\otimes \H^{T}\W^{T})   -  \C^{(n\times r)} (\I_{r}\otimes \X_{\textup{data}})^{T} \right],\\
\nabla_{\vect(\Beta)}\nabla_{\vect(\H)^{T}} \, L(\bZ) &= \nabla_{\vect(\bGamma)}\nabla_{\vect(\H)^{T}} \, L(\bZ) =  O,\\
\nabla_{\vect(\bGamma)}\nabla_{\vect(\Beta)^{T}} \, L(\bZ) &= (\I_{r}\otimes \X_{\textup{aux}}) \diag\left( \ddot{\H}(y_{1},\a_{1}),\dots,\ddot{\H}(y_{n},\a_{n})\right)  \C^{(n,\kappa)} (\I_{\kappa}\otimes \W^{T}\X_{\textup{data}})^{T}.
\end{align}
\end{lemma}

\begin{proof}
Setting $\nu=0$ in \eqref{eq:f_SDL_pf0}, we have $L(\bZ)=f(\A,\B,\bGamma)$. Recall the gradients of $f(\A,\B,\bGamma)$ in \eqref{eq:SDL_filt_gradients}. By using the chain rule and noting that $\A[:,j]=\W \Beta[:,j]$, we can compute
\begin{align}
\nabla_{\W} f(\A,\B,\bGamma) &= \left( \sum_{s=1}^{n} \sum_{j=1}^{\kappa} \frac{\partial f(\A,\B,\bGamma)}{\partial \A[:,j]} \frac{\partial \A[:,j]}{\partial \W}\right) + 2\xi(\W\H-\X_{\textup{data}})\H^{T}  \\
&= \sum_{s=1}^{n}\sum_{j=1}^{\kappa} \x_{s} \dot{h}_{j}(y_{s},\a_{s}) \Beta[:,j]^{T} + 2\xi(\W\H-\X_{\textup{data}})\H^{T} \\
&= \X_{\textup{data}} \K^{T} \Beta^{T} +  2\xi(\W\H-\X_{\textup{data}})\H^{T},
\end{align}
By a similar computation, we can also compute, for each $j=1,\dots,\kappa$,
\begin{align}
\nabla_{\Beta[:,j]} L(\bZ) &= \sum_{s=1}^{n}  \frac{\partial \A[:,j]}{\partial \Beta[:,j]}\frac{\partial f(\A,\B,\bGamma)}{\partial \A[:,j]}  = \sum_{s=1}^{n} \W^{T}  \dot{h}_{j}(y_{s},\a_{s}) \x_{s},
\end{align}
so we get $\nabla_{\Beta} f(\A,\B,\bGamma)  = \W^{T} \X_{\textup{data}}  \K^{T}$. A similar computation shows the remaining two gradients.

Next, recall the relations for vectorizing product of matrices: for $A\in \R^{a\times b}$, $B\in \R^{b\times c}$, and $C\in \R^{c\times d}$,
\begin{align}
\vect(AB) &= (\I_{c}\otimes A) \vect(B) = (B^{T} \otimes \I_{a}) \vect(A), \\
\vect(ABC) &= (C^{T}\otimes A)\vect(B) = (\I_{d}\otimes AB) \vect(C) = (C^{T}B^{T}\otimes \I_{a}) \vect(A).
\end{align}
From this the previous calculation yields
\begin{align}
\nabla_{\vect(\W)} f(\A,\B,\bGamma)
&= \vect\left( \X_{\textup{data}} \K^{T} \Beta^{T} \right) +  2\xi \vect(\W\H\H^{T}) - 2\xi \vect(\X_{\textup{data}} \H^{T}) \\
&= (\Beta \otimes \X_{\textup{data}}) \vect(\K^{T}) + 2\xi (\H\H^{T} \otimes \I_{p} ) \vect(\W)  - 2\xi \vect(\X_{\textup{data}} \H^{T}).
\end{align}
Note that $\vect(\K^{T})^{T}=(\C^{(\kappa,n)} \vect(\K))^{T}=\vect(\K)^{T} \C^{(n,\kappa)}$. Hence we get
\begin{align}
\nabla_{\vect(\W)}\nabla_{\vect(\W)^{T}} f(\A,\B,\bGamma)
&= \nabla_{\vect(\W)}\left( \vect(\K)^{T} \C^{(n,\kappa)} (\Beta \otimes \X_{\textup{data}})^{T}  + 2\xi  \vect(\W)^{T}  (\H\H^{T} \otimes \I_{p} ) - 2\xi \vect(\X_{\textup{data}} \H^{T})^{T} \right)\\
&= (\Beta \otimes \X_{\textup{data}}) \diag\left( \ddot{\H}(y_{1},\a_{1}),\dots,\ddot{\H}(y_{n},\a_{n})\right) \C^{(n,\kappa)} (\Beta \otimes \X_{\textup{data}})^{T} +2\xi (\H\H^{T} \otimes \I_{p} ).
\end{align}
Similarly, we can compute
\begin{align}
\nabla_{\vect(\Beta)}\nabla_{\vect(\Beta)^{T}} f(\A,\B,\bGamma) &=\nabla_{\vect(\Beta)} \vect(\W^{T} \X_{\textup{data}}  \K^{T})^{T}\\
&= \nabla_{\vect(\Beta)} \vect(\K)^{T} (\I_{\kappa}\otimes \W^{T}\X_{\textup{data}})^{T} \\
&=(\I_{\kappa}\otimes \W^{T}\X_{\textup{data}}) \diag\left( \ddot{\H}(y_{1},\a_{1}),\dots,\ddot{\H}(y_{n},\a_{n})\right)   (\I_{\kappa}\otimes \W^{T}\X_{\textup{data}})^{T}.
\end{align}
Also note that
\begin{align}
\nabla_{\vect(\bGamma)}\nabla_{\vect(\bGamma)^{T}} f(\A,\B,\bGamma) &= \nabla_{\vect(\bGamma)} \vect(\X_{\textup{aux}} \K^{T})^{T} \\
&= \nabla_{\vect(\bGamma)} \vect(\K)^{T} (\I_{r}\otimes \X_{\textup{aux}})^{T} \\
&= (\I_{r}\otimes \X_{\textup{aux}}) \diag\left( \ddot{\H}(y_{1},\a_{1}),\dots,\ddot{\H}(y_{n},\a_{n})\right)  (\I_{r}\otimes \X_{\textup{aux}})^{T}.
\end{align}
Similarly, we get
\begin{align}
\nabla_{\vect(\H)}\nabla_{\vect(\H)^{T}} f(\A,\B,\bGamma) &= \nabla_{\vect(\H)} \left( 2\xi \vect(\W^{T}\W\H)^{T} - 2\xi \vect(\W^{T}\X_{\textup{data}})^{T}\right) \\
&= 2\xi \nabla_{\vect(\H)} \vect(\H)^{T} (\I_{n}\otimes \W^{T}\W)  = (\I_{n}\otimes \W^{T}\W).
\end{align}

Next, we compute the off-diagonal terms in the Hessian of $f$. First, we compute
\begin{align}
\nabla_{\vect(\Beta)}\nabla_{\vect(\W)^{T}} f(\A,\B,\bGamma) &= \nabla_{\vect(\Beta)} \vect(\X_{\textup{data}} \K^{T} \Beta^{T})^{T} \\
&= \left( \frac{\partial}{\partial \vect(\Beta) } \vect(\Beta)^{T}\C^{(\kappa, r)} \right) (\I_{r}\otimes \X_{\textup{data}} \K^{T} )^{T} + \left( \frac{\partial}{\partial \vect(\Beta)}  \C^{(\kappa\times n)} \vect(\K^{T}) \right)  (\Beta \otimes \X_{\textup{data}})^{T}   \\
&= \C^{(\kappa, r)} (\I_{r}\otimes \X_{\textup{data}} \K^{T} )^{T} + \\
&\qquad + \C^{(\kappa, n)} (\I_{\kappa}\otimes \W^{T}\X_{\textup{data}})\diag\left( \ddot{\H}(y_{1},\a_{1}),\dots,\ddot{\H}(y_{n},\a_{n})\right)  (\Beta \otimes \X_{\textup{data}})^{T}.
\end{align}
Second, note that $\nabla_{\vect(\bGamma)}\nabla_{\vect(\W)^{T}} f(\A,\B,\bGamma)=O$. Third, for the forthcoming computation, we claim that
\begin{align}
\nabla_{\vect(\H)} \vect(\H\H^{T})^{T} = (\H^{T}\otimes \I_{r}) + (\I_{r}\otimes \H^{T}).
\end{align}
One can directly verify the above when $\H$ consists of a single column, and the general case can be easily obtained from there. Also note that using the commutation matrix, we can write $\vect(\H^{T})^{T}=(\C^{(r,n)} \vect(\H))^{T}=\vect(\H)^{T} \C^{(n,r)}$. Now observe that
\begin{align}
\nabla_{\vect(\H)}\nabla_{\vect(\W)^{T}} f(\A,\B,\bGamma) &= 2\xi \nabla_{\vect(\H)} \left[ \vect(\W\H\H^{T}) - \vect(\X_{\textup{data}} \H^{T}) \right]^{T} \\
&= 2\xi \nabla_{\vect(\H)} \left[ \vect(\H\H^{T})^{T} (\I_{r}\otimes \W)^{T}  -  \vect(\H^{T})^{T}(\I_{r}\otimes \X_{\textup{data}})^{T}  \right] \\
&=   2\xi \left( \nabla_{\vect(\H)}  \vect(\H\H^{T})^{T} \right) (\I_{r}\otimes \W)^{T}  -  \left( \nabla_{\vect(\H)}  \vect(\H^{T})^{T} \right) (\I_{r}\otimes \X_{\textup{data}})^{T} \\
&=  2\xi \left[ \left( (\H^{T}\otimes \I_{r}) + (\I_{r}\otimes \H^{T}) \right) (\I_{r}\otimes \W)^{T}  -  \C^{(n\times r)} (\I_{r}\otimes \X_{\textup{data}})^{T} \right].
\end{align}
Then one can use the mixed-product property to further simplify the last expression as in the assertion. Fourth, noting that $\vect(\K^{T})^{T}=(\C^{(\kappa,n)} \vect(\K))^{T}=\vect(\K)^{T} \C^{(n,\kappa)}$, we can compute
\begin{align}
\nabla_{\vect(\bGamma)}\nabla_{\vect(\Beta)^{T}} f(\A,\B,\bGamma) &= \nabla_{\vect(\bGamma)} \vect(\W^{T} \X_{\textup{data}}  \K^{T})^{T} \\
&= \nabla_{\vect(\bGamma)} \vect(\K^{T})^{T} (\I_{\kappa}\otimes \W^{T}\X_{\textup{data}})^{T}\\
&= \nabla_{\vect(\bGamma)} \vect(\K)  \C^{(n,\kappa)} (\I_{\kappa}\otimes \W^{T}\X_{\textup{data}})^{T} \\
&= (\I_{r}\otimes \X_{\textup{aux}}) \diag\left( \ddot{\H}(y_{1},\a_{1}),\dots,\ddot{\H}(y_{n},\a_{n})\right)  \C^{(n,\kappa)} (\I_{\kappa}\otimes \W^{T}\X_{\textup{data}})^{T}.
\end{align}
The remaining zero-second derivatives are easy to see.
\end{proof}

\begin{remark}[Derivatives for the feature-based SDL objective in separate variables]
\label{rmk:feat_SDL_gradient}
Arguing similarly as in the proof of Lemma \ref{lem:SDL_filt_BCD_derivatives}, we can compute the derivatives of the feature-based SDL objective in separate variables as follows. Let $L(\bZ)$ denote the objective of the feature-based SDL in \eqref{eq:ASDL_1}. Suppose \ref{assumption:A4} holds. Recall $\dot{\h}$ defined in \eqref{eq:Hddot_def}. Let $\a_{s}:=\Beta^{T}\h_{s}+\bGamma^{T}\x'_{s}$ for $s=1,\dots,n$, where $\H=[\h_{1},\dots,\h_{n}]\in \R^{r\times n}$ being the code matrix. Let $\K:=[\dot{\h}(y_{1},\a_{1}),\dots, \dot{\h}(y_{1},\a_{1})]\in \R^{\kappa\times n}$. Then we
\begin{align}\label{eq:SDL_feat_BCD_derivatives1}
\nabla_{\W} \,  L(\bZ) &= 2\xi(\W\H-\X_{\textup{data}})\H^{T},\qquad
\nabla_{\Beta} \, L(\bZ)  = \H  \K^{T} \\
\nabla_{\bGamma} \, L(\bZ) &=  \X_{\textup{aux}}\K^{T}, \qquad  \nabla_{\H} \, L(\bZ) =  \Beta \K +   2\xi \W^{T}(\W\H-\X_{\textup{data}}).
\end{align}
\end{remark}

\begin{lemma}\label{lem:SDL_BCD_smooth}
\label{lem:BCD_hypothesis}
Assume the hypothesis of Theorem \ref{thm:SDL_BCD} is true. Then the loss function $L$ in \eqref{eq:ASDL_1} is convex in each block coordinates $\W$, $\H$, $\Beta$, and $\bGamma$. Furthermore, its gradient is continuous and $M$-Lipschitz for some $M>0$ on the admissible parameter space.
\end{lemma}

\begin{proof}
Notice that the diagonal terms in the Hessian given in  Lemma \ref{lem:SDL_filt_BCD_derivatives} are positive semidefinite. (An explicit lower bound on the eigenvalues can also be computed.) This is enough to conclude that $L$ is convex in each factor $\W$, $\H$, $\Beta$, and $\bGamma$ while all the other three are held fixed (i.e., $L$ is multiconvex).  The second part of the assertion follows easily from the first derivative computations in Lemma \ref{lem:SDL_filt_BCD_derivatives} and the compactness assumption in Theorem \ref{thm:SDL_BCD}. %The second part of the assertion also follows easily from Lemma \ref{lem:SDL_filt_BCD_derivatives}, \ref{assumption:A5_cpt}, and the dominated convergence theorem.
\end{proof}

Now we are ready to derive Theorem \ref{thm:SDL_BCD}.

\begin{proof}[\textbf{Proof of Theorem \ref{thm:SDL_BCD}}]
The result would immediately follow the main result in [44]. In order to apply the result, we need to verify that 1) the filter-based SDL loss function $L$ in \eqref{eq:ASDL_1} is multiconvex and 2) the gradient of $L$ is $M$-Lipschitz for some constant $M>0$. Under the assumption, \ref{assumption:A4} and the assumed compactness of the parameter space, both of these hypotheses are verified in Lemma \ref{lem:BCD_hypothesis}. This shows the assertion.
\end{proof}

\section{Proof of Theorem \ref{thm:SDL_BCD_STAT_filt}}
\label{section:appendix_SDL_BCD_STAT}

Throughout this section, let $\mathcal{L}(\bZ)=\mathcal{L}_{n}(\W,\h,\Beta,\bGamma)$ denote the objective in \eqref{eq:SDL_likelihood_BCD_filter}. Denote
\begin{align}\label{eq:def_L_bar_STAT}
\bar{\mathcal{L}}(\bZ) := \E_{(\x,\x',y)}\left[ \mathcal{L}_{1}(\bZ) \right].
\end{align}

\begin{lemma}[Derivatives of the filter-based SDL objective in separate variables]
\label{lem:SDL_filt_BCD_STAT_derivatives}
Let $\mathcal{L}(\bZ)=\mathcal{L}_{n}(\W,\h,\Beta,\bGamma)$ denote the objective in \eqref{eq:SDL_likelihood_BCD_filter}. Suppose \ref{assumption:A4} holds. Recall $\dot{\h}$ and $\ddot{\H}$ defined in \eqref{eq:Hddot_def}. Let $\a_{s}:=\Beta^{T}\W^{T}\x_{s}+\bGamma^{T}\x'_{s}$ for $s=1,\dots,n$ and $\K:=[\dot{\h}(y_{1},\a_{1}),\dots, \dot{\h}(y_{1},\a_{1})]\in \R^{\kappa\times n}$. Also denote $\H=[\h,\dots,\h]\in \R^{r\times n}$. Then we have
\begin{align}
\nabla_{\W} \,  \mathcal{L}(\bZ) &= \X_{\textup{data}} \K^{T} \Beta^{T} +  2\xi(\W\H-\X_{\textup{data}})\H^{T} + 2\nu\W,\qquad
\nabla_{\Beta} \, \mathcal{L}(\bZ)  = \W^{T} \X_{\textup{data}}  \K^{T} + 2\nu \Beta, \\
\nabla_{\bGamma} \, \mathcal{L}(\bZ) &=  \X_{\textup{aux}}\K^{T} + 2\nu \bGamma, \qquad  \nabla_{\h} \, \mathcal{L}(\bZ) =  2\xi \W^{T}(n\W\h-\sum_{s=1}^{n}\x_{s}) + 2\nu \h.
\end{align}
Furthermore, for diagonal terms in the Hessian, we have
\begin{align}
\nabla_{\vect(\W)}\nabla_{\vect(\W)^{T}} \, \mathcal{L}(\bZ)
&= (\Beta \otimes \X_{\textup{data}}) \diag\left( \ddot{\H}(y_{1},\a_{1}),\dots,\ddot{\H}(y_{n},\a_{n})\right) \C^{(n,\kappa)} (\Beta \otimes \X_{\textup{data}})^{T} +2\xi (\H\H^{T} \otimes \I_{p} ) + 2\nu \I_{pr}, \\
\nabla_{\h}\nabla_{\h^{T}} \, \mathcal{L}(\bZ) & = 2\xi n\W^{T}\W + 2\nu \I_{r}, \\
\nabla_{\vect(\Beta)}\nabla_{\vect(\Beta)^{T}} \, \mathcal{L}(\bZ)
&=(\I_{\kappa}\otimes \W^{T}\X_{\textup{data}}) \diag\left( \ddot{\H}(y_{1},\a_{1}),\dots,\ddot{\H}(y_{n},\a_{n})\right)   (\I_{\kappa}\otimes \W^{T}\X_{\textup{data}})^{T} + 2\nu \I_{r\kappa}, \\
\nabla_{\vect(\bGamma)}\nabla_{\vect(\bGamma)^{T}} \, \mathcal{L}(\bZ)
&= (\I_{r}\otimes \X_{\textup{aux}}) \diag\left( \ddot{\H}(y_{1},\a_{1}),\dots,\ddot{\H}(y_{n},\a_{n})\right)  (\I_{r}\otimes \X_{\textup{aux}})^{T} +  2\nu \I_{q\kappa}
\end{align}
Lastly, for the off-diagonal terms in the Hessian, we have
\begin{align}
\nabla_{\vect(\Beta)}\nabla_{\vect(\W)^{T}} \, \mathcal{L}(\bZ)
&=  (\I_{\kappa}\otimes \W^{T}\X_{\textup{data}})\diag\left( \ddot{\H}(y_{1},\a_{1}),\dots,\ddot{\H}(y_{n},\a_{n})\right)  (\Beta \otimes \X_{\textup{data}})^{T}  +  \C^{(\kappa, r)} (\I_{\r}\otimes \X_{\textup{data}} \K^{T} )^{T}, \\
\nabla_{\vect(\bGamma)}\nabla_{\vect(\W)^{T}} \, \mathcal{L}(\bZ) &= O,\\
\nabla_{\h}\nabla_{\vect(\W)^{T}} \, \mathcal{L}(\bZ) &= 2\xi \left[ n(\h^{T}\otimes \W^{T}) + n(\I_{r}\otimes \h^{T}\W^{T})   -  (\mathbf{1}_{1\times n}\otimes \I_{r}) \C^{(n\times r)} (\I_{r}\otimes \X_{\textup{data}})^{T} \right],\\
\nabla_{\vect(\Beta)}\nabla_{\h^{T}} \, \mathcal{L}(\bZ) &= \nabla_{\vect(\bGamma)}\nabla_{\h^{T}} \, \mathcal{L}(\bZ) =  O,\\
\nabla_{\vect(\bGamma)}\nabla_{\vect(\Beta)^{T}} \, \mathcal{L}(\bZ) &= (\I_{r}\otimes \X_{\textup{aux}}) \diag\left( \ddot{\H}(y_{1},\a_{1}),\dots,\ddot{\H}(y_{n},\a_{n})\right)  \C^{(n,\kappa)} (\I_{\kappa}\otimes \W^{T}\X_{\textup{data}})^{T}.
\end{align}
\end{lemma}

\begin{proof}
For $\H=[\h,\dots,\h]$, note that
\begin{align}
\nabla_{\h} \vect(\H)^{T} = \mathbf{1}_{1\times n}\otimes \I_{r} ,\qquad    \nabla_{\h} \vect(\H\H^{T})^{T} = n(\h^{T}\otimes \I_{r}) + n(\I_{r}\otimes \h^{T}).
\end{align}
Then the assertion follows from similar computations as in the proof of Lemma \ref{lem:SDL_filt_BCD_derivatives}.
\end{proof}

\begin{lemma}\label{lem:SDL_BCD_STAT_smooth}
\label{lem:BCD_STAT_hypothesis}
Let $\mathcal{L}_{n}(\bZ)=\mathcal{L}_{n}(\W,\h,\Beta,\bGamma)$ denote the objective in \eqref{eq:SDL_likelihood_BCD_filter}. Assume the hypothesis of Theorem \ref{thm:SDL_BCD_STAT_filt} holds. Then $\mathcal{L}_{n}$ is convex in each block coordinates $\W$, $\h$, $\Beta$, and $\bGamma$ for all $\nu\ge 0$. Furthermore, its gradient is continuous and $M$-Lipschitz for some $M>0$ on the admissible parameter space.
\end{lemma}

\begin{proof}
The argument is identical to the proof Lemma \ref{lem:BCD_hypothesis} using Lemma \ref{lem:SDL_filt_BCD_STAT_derivatives} instead of Lemma \ref{lem:SDL_filt_BCD_derivatives}. For the multi-convexity part, recall that $\ddot{\H}(y_{s},\a_{s})$'s are positive definite due to \ref{assumption:A4} and the commutation matrices $\C^{(a,b)}$ are also positive definite.
\end{proof}

\begin{lemma}[Derivatives of the expected filter-based SDL objective]
\label{lem:SDL_filt_BCD_STAT_derivatives}
Suppose a single data $(\x,\x',y)$ is sampled according to the generative model \eqref{eq:SDL_prob_BCD_filt} and let $\bar{\mathcal{L}}$ be as in \eqref{eq:def_L_bar_STAT}. Assume the hypothesis of Theorem \ref{thm:SDL_BCD_STAT_filt} holds.
Recall $\dot{\h}$ and $\ddot{\H}$ defined in \eqref{eq:Hddot_def}. Let $\a:=\Beta^{T}\W^{T}\x+\bGamma^{T}\x'$ and $\xi:=(2\sigma^{2})^{-1}$. Then we have
\begin{align}
\nabla_{\W} \,  \bar{\mathcal{L}}(\bZ) &= \E[\x \, \dot{\h}(y,\a)^{T}] \Beta^{T} +  2\xi(\W\h-\W^{\star}\h^{\star})\h^{T} + 2\nu \W,\qquad
\nabla_{\Beta} \, \bar{\mathcal{L}}(\bZ)  = \W^{T} \E[\x  \,\dot{\h}(y,\a)^{T}] + 2\nu \Beta\\
\nabla_{\bGamma} \, \bar{\mathcal{L}}(\bZ) &=  \E[\x' \dot{\h}(y,\a)^{T}] + 2\nu \bGamma \qquad  \nabla_{\h} \, \bar{\mathcal{L}}(\bZ) =  2\xi \W^{T}(\W\h-\W^{\star}\h^{\star}) + 2\nu \h.
\end{align}
Furthermore, for diagonal terms in the Hessian, we have
\begin{align}
\nabla_{\vect(\W)}\nabla_{\vect(\W)^{T}} \, \bar{\mathcal{L}}(\bZ)
&=\E\left[ (\Beta \otimes \x)  \ddot{\H}(y,\a) (\Beta \otimes \x)^{T} \right]+2\xi (\h\h^{T} \otimes \I_{p} ) + 2\nu \I_{pr} \, \\
\nabla_{\h}\nabla_{\h^{T}} \, \bar{\mathcal{L}}(\bZ) & = 2\xi \W^{T}\W +  2\nu \I_{r},\\
\nabla_{\vect(\Beta)}\nabla_{\vect(\Beta)^{T}} \, \bar{\mathcal{L}}(\bZ)
&=\E\left[(\I_{\kappa}\otimes \W^{T}\x)  \ddot{\H}(y,\a)  (\I_{\kappa}\otimes \W^{T}\x)^{T}\right] + 2\nu \I_{r\kappa}, \\
\nabla_{\vect(\bGamma)}\nabla_{\vect(\bGamma)^{T}} \, \bar{\mathcal{L}}(\bZ)
&= \E\left[(\I_{r}\otimes \x') \ddot{\H}(y,\a) (\I_{r}\otimes \x')^{T}\right]+ 2\nu \I_{q\kappa}.
\end{align}
Lastly, for the off-diagonal terms in the Hessian, we have
\begin{align}
\nabla_{\vect(\Beta)}\nabla_{\vect(\W)^{T}} \, \bar{\mathcal{L}}(\bZ)
&= \E\left[  (\I_{\kappa}\otimes \W^{T}\x)  \ddot{\H}(y,\a) (\Beta \otimes \x)^{T}  +  \C^{(\kappa\times r)} (\I_{r}\otimes \x \, \dot{\h}(y, \a)^{T} )^{T} \right], \\
\nabla_{\vect(\bGamma)}\nabla_{\vect(\W)^{T}} \, \bar{\mathcal{L}}(\bZ) &= O,\\
\nabla_{\h}\nabla_{\vect(\W)^{T}} \, \bar{\mathcal{L}}(\bZ) &= 2\xi \left[ (\h^{T}\otimes \W^{T}) + (\I_{r}\otimes (\W\h - \W^{\star}\h^{\star})^{T})  \right],\\
\nabla_{\vect(\Beta)}\nabla_{\h^{T}} \, \bar{\mathcal{L}}(\bZ) &= \nabla_{\vect(\bGamma)}\nabla_{\h^{T}} \, \bar{\mathcal{L}}(\bZ) =  O,\\
\nabla_{\vect(\bGamma)}\nabla_{\vect(\Beta)^{T}} \, \bar{\mathcal{L}}(\bZ) &= \E\left[ (\I_{r}\otimes \x')  \ddot{\H}(y,\a) (\I_{\kappa}\otimes \W^{T}\x)^{T} \right].
\end{align}
\end{lemma}

\begin{proof}
According to Lemmas \ref{lem:SDL_filt_BCD_STAT_derivatives} and \ref{lem:SDL_BCD_STAT_smooth}, $\mathcal{L}$ is twice continuously differentiable and both $\nabla \mathcal{L}$ and $\nabla^{2} \mathcal{L}$ are bounded within the (compact) parameter space. Hence by the monotone convergence theorem,
\begin{align}
\nabla \bar{\mathcal{L}} = \E\left[ \nabla \mathcal{L}\right], \qquad \nabla^{2} \bar{\mathcal{L}} = \E\left[ \nabla^{2} \mathcal{L}\right].
\end{align}
Hence we can simply specialize the derivatives of $\mathcal{L}$ we computed in Lemma \ref{lem:SDL_filt_BCD_STAT_derivatives} for the single sample case $n=1$ and then take the expectation. In doing so, we use the fact that $\C^{(1\times a)}=\I_{a}$, where $\C^{(a,b)}$ denotes the commutation matrix defined above the statement of Lemma \ref{lem:SDL_filt_BCD_derivatives}. %The derivation is straightforward and details are omitted.
\end{proof}

\begin{lemma}\label{lem:SDL_BCD_avg_smooth}
\label{lem:BCD_STAT_hypothesis}
Suppose a single data $(\x,\x',y)$ is sampled according to the generative model \eqref{eq:SDL_prob_BCD_filt} with true parameters $\bZ^{\star}=[\W^{\star}, \H^{\star}, \Beta^{\star}, \bGamma^{\star}]$ and $\blambda^{\star}$. Let $\bar{\mathcal{L}}$ be as in \eqref{eq:def_L_bar_STAT} and assume the hypothesis of Theorem \ref{thm:SDL_BCD_STAT_filt} holds. Then the following hold:
\begin{description}
\item[(i)] $\bar{\mathcal{L}}$ is convex in each block coordinates $\W$, $\h$, $\Beta$, and $\bGamma$ for all $\nu\ge 0$.

\item[(ii)] $\nabla \bar{\mathcal{L}}$ continuous and $M$-Lipschitz for some $M>0$ on the admissible parameter space.

\item[(iii)] $\nabla^{2} \bar{\mathcal{L}}(\bZ^{\star})$ is positive definite if $\nu > \lambda_{+}$, where
\begin{align}
\lambda_{+}:= \frac{1}{2}\max\left( \begin{matrix} 2\xi \lVert \h^{\star}  \rVert_{2} \, \lVert  \W^{\star} \rVert_{2} + \alpha^{+} \lVert \Beta^{\star} \rVert_{2}  \lVert \W^{\star} \rVert_{2} \E\left[  \Vert \x\x^{T} \rVert_{2} \right] +  \gamma_{\max}\sigma \sqrt{2p\pi}  \\
\hspace{3cm} - \alpha^{-} \, \lambda_{\min}(\Beta^{\star} (\Beta^{\star})^{T} ) \E\left[   \lambda_{\min}( \x\x^{T}  ) \right] - 2\xi \lambda_{\min} (\h^{\star}(\h^{\star})^{T} ), \\[10pt]
2\xi \lVert \h^{\star}  \rVert_{2} \, \lVert  \W^{\star} \rVert_{2} - 2\xi (\W^{\star})^{T}\W^{\star} - 2\xi \lambda_{\min}( (\W^{\star})^{T}\W^{\star}), \\[10pt]
\alpha^{+} \lVert \Beta^{\star} \rVert_{2}  \lVert \W^{\star} \rVert_{2} \E\left[  \Vert \x\x^{T} \rVert_{2} \right] +  \gamma_{\max}\sigma \sqrt{2p\pi}  +\alpha^{+} \E\left[ \lVert \x'(\W^{\star})^{T}\x \rVert_{2}\right] \\
\hspace{4cm} - \alpha^{-} \, \E\left[  \lambda_{\min}( (\W^{\star})^{T} \x\x^{T} \W^{\star} ) \right], \\[10pt]
\alpha^{+} \E\left[ \lVert \x'(\W^{\star})^{T}\x \rVert_{2}\right] - \lambda_{\min}(\x'(\x')^{T})
\end{matrix}\right) \label{eq:nu_min_BCD_STAT_explicit}.
\end{align}
\end{description}
\end{lemma}

\begin{proof}
Parts \textbf{(i)} and \textbf{(ii)} follow easily from Lemma \ref{lem:SDL_filt_BCD_STAT_derivatives} as in the proof of Lemma \ref{lem:BCD_STAT_hypothesis}. Now we argue for \textbf{(iii)}. Denote $\xi=(2\sigma^{2})^{-1}$ and  $\a^{\star}=(\W^{\star}\Beta^{\star})^{T}\x+\bGamma^{\star}\x' $. Recall that according to Lemma \ref{lem:SDL_filt_BCD_STAT_derivatives}, we can write the Hessian of $\bar{\mathcal{L}}$ at the true parameter $\bZ^{\star}$ as the $4\times 4$ block matrix $(A_{ij})_{1\le i,j\le 4} + 2\nu \I$ in \eqref{eq:expected_hessian_block}. Recall \ref{assumption:A4}. The diagonal blocks are given by
\begin{align}
&=\E\left[ (\Beta^{\star} \otimes \x)  \ddot{\H}(y,\a^{\star}) (\Beta^{\star} \otimes \x)^{T} \right]+2\xi (\h^{\star}(\h^{\star})^{T} \otimes \I_{p} ) \\
&\qquad \succeq \left( \alpha^{-} \, \lambda_{\min}(\Beta^{\star} (\Beta^{\star})^{T} ) \E\left[   \lambda_{\min}( \x\x^{T}  ) \right] + 2\xi \lambda_{\min} (\h^{\star})(\h^{\star})^{T}  \right) \I_{pr}   \\
A_{22} & = 2\xi (\W^{\star})^{T}\W^{\star} \succeq 2\xi \lambda_{\min}( (\W^{\star})^{T}\W^{\star}) \I_{r},\\
A_{33}&=\E\left[(\I_{\kappa}\otimes (\W^{\star})^{T}\x)  \ddot{\H}(y,\a^{\star})  (\I_{\kappa}\otimes (\W^{\star})^{T}\x)^{T}\right] \succeq \alpha^{-} \E\left[ \lambda_{\min}( (\W^{\star})^{T}\x \x^{T} \W^{\star} ) \right] \I_{r\kappa}, \\
&= \E\left[(\I_{r}\otimes \x') \ddot{\H}(y,\a^{\star}) (\I_{r}\otimes \x')^{T}\right] \succeq \alpha^{-} \lambda_{\min}(\x'(\x')^{T}) \I_{rq},
\end{align}
where the off-diagonal blocks are given by
\begin{align}
A_{21} &=   2\xi  (\h^{\star})^{T}\otimes (\W^{\star})^{T} \\
A_{31} &=  \E\left[  (\I_{\kappa}\otimes \W^{T}\x)  \ddot{\H}(y,\a^{\star}) (\Beta \otimes \x)^{T}  +  \C^{(\kappa\times r)} (\I_{r}\otimes \x \, \dot{\h}(y,\a^{\star})^{T} )^{T} \right]  \\
A_{43} &=  \E\left[ (\I_{r}\otimes \x')  \ddot{\H}(y,\a^{\star}) (\I_{\kappa}\otimes \W^{T}\x)^{T} \right].
\end{align}
If $\nu$ is large enough so that the following condition is satisfied
\begin{align}
\lambda_{\min}(A_{ii}) + 2\nu > \sum_{j\ne i} \lVert A_{ij} \rVert_{2} \quad \forall 1\le i \le 4,
\end{align}
then the Hessian $\nabla^{2} \bar{\mathcal{L}}$ is block diagonally dominant and is positive definite (see [22]). Thus it suffices to take
\begin{align}\label{eq:lambda_lower_bd0}
\nu &> \frac{1}{2}\max\left( \begin{matrix} \lVert A_{12} \rVert_{2} +\lVert A_{13} \rVert_{2}   - \lambda_{\min}(A_{11}), \\
\lVert A_{12} \rVert_{2}- \lambda_{\min}(A_{22}), \\
\lVert A_{13} \rVert_{2} +\lVert A_{34} \rVert_{2} - \lambda_{\min}(A_{33}), \\
\lVert A_{34} \rVert_{2} - \lambda_{\min}(A_{44})
\end{matrix}\right).
\end{align}
Note that (see \ref{assumption:A4} for the definition of $\gamma_{\max}$ and $\alpha^{\pm}$)
\begin{align}
\lVert \E\left[  \x \, \dot{\h}(y,\a^{\star})^{T} \right] \rVert_{2} &\le \lVert \E\left[  \W^{\star} \h^{\star} \dot{\h}(y,\a^{\star})^{T} \right] \rVert_{2} + \lVert \E\left[ \beps \dot{\h}(y,\a^{\star})^{T} \right] \rVert_{2} \\
&\le  \lVert \W^{\star} \h^{\star} \E\left[  \dot{\h}(y,\a^{\star})^{T} \right] \rVert_{2} + \lVert \E\left[ \beps \dot{\h}(y,\a^{\star})^{T} \right] \rVert_{2} \\
&= \lVert \E\left[ \beps \,\dot{\h}(y, \a^{\star})^{T} \right] \rVert_{2} \\
&\le \gamma_{\max}\sigma \sqrt{2p\pi}.
\end{align}
Using this, we get
\begin{align}
\lVert A_{12} \rVert_{2} &= \lVert A_{21} \rVert_{2} =2\xi \lVert \h^{\star}  \rVert_{2} \, \lVert  \W^{\star} \rVert_{2}, \\
\lVert A_{13} \rVert_{2} &\le \alpha^{+} \lVert \Beta^{\star} \rVert_{2}  \lVert \W^{\star} \rVert_{2} \E\left[  \Vert \x\x^{T} \rVert_{2} \right] +  \gamma_{\max}\sigma \sqrt{2p\pi}, \\
\lVert A_{43} \rVert_{2} &\le \alpha^{+} \E\left[ \lVert \x'(\W^{\star})^{T}\x \rVert_{2}\right].
\end{align}
Using these upper bounds on the operator norm of the off-diagonal blocks and the lower bounds on the eigenvalues of the diagonal blocks above, the lower bound in \eqref{eq:lambda_lower_bd0} can be lower bounded by $\lambda_{+}$ defined in the assertion.
\end{proof}

\section{Proof of Theorems \ref{thm:SDL_LPGD_STAT} and \ref{thm:SDL_LPGD_STAT_feat}}

We first recall the following standard concentration bounds:

\begin{lemma}[Generalized Hoeffding's inequality for sub-gaussian variables]
\label{lem:hoeffding}
Let $X_{1},\dots,X_{n}$ denote i.i.d. random vectors in $\R^{d}$ such that $\E[X_{k}[i]^{2}/K^{2}]\le 2$ for some constant $K>0$ for all $1\le k \le n$ and $1\le i \le d$. Fix a vector $\a=(a_{1},\dots,a_{n})^{T}\in \R^{n}$. Then for each $t>0$,
\begin{align}
\P\left( \left\lVert \sum_{k=1}^{n} a_{k}X_{k} \right\rVert_{1} > t \right) \le 2d \exp\left( \frac{-t^{2}}{K^{2} d^{2}\lVert \a  \rVert_{2}^{2}} \right)
\end{align}
\end{lemma}

\begin{proof}
Follows from \cite[Thm 2.6.2]{vershynin2018high} and using a union bound over $d$ coordinates.
\end{proof}

\begin{lemma}(2-norm of matrices with independent sub-gaussian entries)
\label{lem:matrix_norm_bd}
Let $\A$ be an $m\times n$ random matrix with independent subgaussian entries $\A_{ij}$ of mean zero. Denote $K$ to be the maximum subgaussian norm of $\A_{ij}$, that is, $K>0$ is the smallest number such that $\E[\exp(\A_{ij})^{2}/K^{2}  ]\le 2$. Then for each $t>0$,
\begin{align}
\P\left( \lVert \A \rVert_{2} \ge 3K(\sqrt{m}+\sqrt{n}+t) \right) \le 2\exp(-t^{2}).
\end{align}
\end{lemma}

\begin{proof}
See \cite[Thm. 4.4.5]{vershynin2018high}
\end{proof}

\begin{proof}[\textbf{Proof of Theorem \ref{thm:SDL_LPGD_STAT}}]
Let $\mathcal{L}_{n}$ denote the $L_{2}$-regularized negative joint negative log likelihood function in \eqref{eq:SDL_likelihood_conv_filter} without the last three terms, and define the  expected loss function $\bar{\mathcal{L}}_{n}(\bZ):= \E_{\beps_{i},\beps_{i}',1\le i \le n}\left[ \mathcal{L}_{n}(\bZ) \right]$. We omit the constant terms in these functions. Define the following gradient mappings of $\bZ^{\star}$ with respect to the empirical $f_{n}$ and the expected $\bar{f}_{n}$ loss functions:
\begin{align}
G(\bZ^{\star}, \tau)=\frac{1}{\tau}\left( \bZ^{\star} - \Pi_{\Param} \left( \bZ^{\star} - \tau \nabla \mathcal{L}_{n}(\bZ^{\star}) \right)\right), \qquad  \bar{G}(\bZ^{\star},\tau):=\frac{1}{\tau}\left( \bZ^{\star} - \Pi_{\Param} \left( \bZ^{\star} - \tau \nabla \bar{\mathcal{L}}_{n}(\bZ^{\star}) \right)\right).
\end{align}
It is elementary to show that the true parameter $\bZ^{\star}$ is a stationary point of $\bar{\mathcal{L}}- \nu (\lVert \A\rVert_{F}^{2} + \lVert \bGamma \rVert_{F}^{2}) $ over $\Param\subseteq \R^{p\times (\kappa+n)}\times \R^{q\times \kappa}$. Hence we have $\bar{G}(\bZ^{\star},\tau)= 2\nu[\A^{\star},O,\bGamma^{\star}]$, so we may write
\begin{align}\label{eq:grad_mapping_compare_stationary}
G(\bZ^{\star}, \tau) &= G(\bZ^{\star}, \tau) - \bar{G}(\bZ^{\star}, \tau) + 2\nu[\A^{\star},O,\bGamma^{\star}] \\
&= \frac{1}{\tau}\left[  \Pi_{\Param}\left( \bZ^{\star}-\tau\nabla \mathcal{L}_{n}(\bZ^{\star}) \right) - \Pi_{\Param}\left( \bZ^{\star}-\tau\nabla \bar{\mathcal{L}}_{n}(\bZ^{\star}) \right) \right] +  2\nu[\A^{\star},O,\bGamma^{\star}]
\end{align}

First, suppose $\bZ^{\star}-\tau \nabla  \mathcal{L}_{n} (\bZ^{\star})\in \Param$ (In particular, this is the case whe $\Param$ equals the whole space). Then we can disregard the projection $\Pi_{\Param}$ in the above display so we get
\begin{align}
G(\bZ^{\star}, \tau) - 2\nu[\A^{\star},O,\bGamma^{\star}] = \nabla \mathcal{L}_{n}(\bZ^{\star}) - \nabla \bar{\mathcal{L}}(\bZ^{\star}) =: [\Delta \X^{\star}, \Delta \bGamma^{\star}].
\end{align}
According to Theorem \ref{thm:SDL_LPGD}, it now suffices show that $G(\bZ^{\star}, \tau)$ above is small with high probability. We use the notation $\U=[\A^{T}, \bGamma^{T}]^{T}$, $\U^{\star}=[(\A^{\star})^{T}, (\bGamma^{\star})^{T}]^{T}$, $\bPhi=[\bphi_{1},\dots,\bphi_{n}]=[\X_{\textup{data}}^{T}, \X_{\textup{aux}}^{T}]^{T}$ (see also the proof of Theorem \ref{thm:SDL_LPGD}). Denote $\a_{s}=\U^{T}\bphi_{s}$ and $\a_{s}^{\star}=(\U^{\star})^{T}\bphi_{s}$ for $s=1,\dots,n$ and introduce the following random quantities
\begin{align}
\mathtt{Q}_{1}:= \sum_{s=1}^{n} \dot{\h}(y_{s},\a_{s}^{\star})\in \R^{\kappa} ,\quad \mathtt{Q}_{2}:=\sum_{s=1}^{n}  \beps_{s}\in \R^{p} ,\quad \mathtt{Q}_{3}:=\sum_{s=1}^{n}  \beps_{s}'\in \R^{q}, \quad \mathtt{Q}_{4}:= [ \beps_{1},\dots,\beps_{n}]  \in \R^{p\times n}.
\end{align}

Recall that
\begin{align}
&\nabla_{\vect(\U)}  \mathcal{L}_{n}(\U,\B) = \left( \sum_{s=1}^{n} \dot{\h}(y_{s},\a_{s}) \otimes \bphi_{s} \right) + 2\nu \vect(\U), \quad \nabla_{\B}  \mathcal{L}_{n}(\U,\B) = \frac{2}{2\sigma^{2}} (\B-\X_{\textup{data}}), \\
&\nabla_{\vect(\U)}  \bar{\mathcal{L}}_{n}(\U,\B) = \left(  \sum_{s=1}^{n} \E\left[ \dot{\h}(y_{s},\a_{s}) \otimes \bphi_{s} \right] \right) + 2\nu \vect(\U), \quad \nabla_{\B}  \bar{\mathcal{L}}_{n}(\U,\B) = \frac{2}{2\sigma^{2}} (\B-\B^{\star}),
\end{align}
where $\dot{\h}$ is defined in \eqref{eq:Hddot_def}. Note that
\begin{align}
\E\left[ \dot{\h}(y_{s},\a_{s}) \,\bigg|\, \bphi_{s} \right]  		&= \left[ \left(\frac{h'(\a[j])}{1+\sum_{c=1}^{\kappa} h(\a[c])} - g_{j}(\a_{s}^{\star})\frac{h'(\a[j])}{h(\a[j])}  \right)_{\a=\a_{s}} \, ; \, j=1,\dots,\kappa \right]\\
&= \left[ \left(\frac{h'(\a[j])}{1+\sum_{c=1}^{\kappa} h(\a[c])} - \frac{h(\a_{s}^{\star}[j])}{1+\sum_{c=1}^{\kappa} h(\a_{s}^{\star}[c])} \frac{h'(\a[j])}{h(\a[j])}  \right)_{\a=\a_{s}} \, ; \, j=1,\dots,\kappa \right],
\end{align}
so the above vanishes when $\a_{s}=\a_{s}^{\star}$. Hence
\begin{align}\label{eq:dot_h_exp_vanish}
\E\left[ \dot{\h}(y_{s},\a_{s}^{\star}) \otimes \bphi_{s} \right] = \E\left[ \E\left[ \dot{\h}(y_{s},\a_{s}^{\star}) \otimes \bphi_{s}\,\bigg|\, \bphi_{s} \right]  \right] =\mathbf{0},
\end{align}
Hence we can compute the following gradients
\begin{align}
\nabla_{\vect(\A)} (\mathcal{L}_{n} - \bar{\mathcal{L}}_{n})(\A,\B,\bGamma)   &=\left(  \sum_{s=1}^{n} \dot{\h}(y_{s},\a_{s}) \otimes \x_{s} \right)  \\
\nabla_{\vect(\bGamma)} (\mathcal{L}_{n} - \bar{\mathcal{L}}_{n})(\A,\B,\bGamma) &=\left( \sum_{s=1}^{n} \dot{\h}(y_{s},\a_{s}) \otimes \x_{s}' \right) \\
\nabla_{\B}  (\mathcal{L}_{n} - \bar{\mathcal{L}}_{n})(\A,\B,\bGamma) &=\frac{2}{2\sigma^{2})} (\B^{\star}-\X_{\textup{data}}) = \frac{2}{2\sigma^{2}} [ \beps_{1},\dots,\beps_{n}] \\
\nabla_{\blambda}  (\mathcal{L}_{n} - \bar{\mathcal{L}}_{n})(\A,\B,\bGamma)  &=\frac{2}{2\sigma^{2}} \sum_{s=1}^{n}\beps_{s}' .
\end{align}
It follows that (recall the definition of $\gamma_{\max}$in \ref{assumption:A4})
\begin{align}
\lVert \nabla_{\A} (\mathcal{L}_{n} - \bar{\mathcal{L}}_{n})(\A^{\star},\B^{\star},\bGamma^{\star})  \rVert_{2}  &= \left\rVert\sum_{s=1}^{n} (\B^{\star}[:,s]+\beps_{s}) \dot{\h}(y_{s},\a_{s}^{\star})^{T} \right\rVert_{2}  \\
& \le \left\rVert  \sum_{s=1}^{n} \B^{\star}[:,s] \dot{\h}(y_{s},\a_{s}^{\star})^{T} \right\rVert_{2}  +\left\rVert \sum_{s=1}^{n} \beps_{s} \dot{\h}(y_{s},\a_{s}^{\star})^{T} \right\rVert_{2}  \\
&\le  \lVert \B^{\star}\rVert_{\infty} \left\rVert  \mathtt{Q}_{1} \right\rVert_{2}  +  \gamma_{\max}  \left\rVert \mathtt{Q}_{2}  \right\rVert_{2}.
\end{align}
Similarly, we have
\begin{align}
\lVert \Delta \bGamma^{\star} \rVert_{F}=	\lVert \nabla_{\bGamma} (\mathcal{L}_{n} - \bar{\mathcal{L}}_{n})(\A^{\star},\B^{\star},\bGamma^{\star})  \rVert_{F} &= \lVert \nabla_{\vect(\bGamma)} (\mathcal{L}_{n} - \bar{\mathcal{L}}_{n})(\A^{\star},\B^{\star},\bGamma^{\star})  \rVert_{2} \\
&\le  q\lVert \blambda^{\star}\rVert_{\infty}  \left\rVert \mathtt{Q}_{1}  \right\rVert_{2}  + q\gamma_{\max}  \left\rVert \mathtt{Q}_{3}  \right\rVert_{2}
\end{align}
Using the fact that $\lVert [A,B] \rVert_{2}\le \lVert A \rVert_{2} +  \lVert B \rVert_{2} $ for two matrices $A,B$ with the same number of rows, we have
\begin{align}\label{eq:SDL_MLE_pf_bd_Q}
\left\rVert \Delta \X^{\star}  \right\rVert_{2} &= \left\lVert  \nabla_{\A} (\mathcal{L}_{n} - \bar{\mathcal{L}}_{n})(\A^{\star},\B^{\star},\bGamma^{\star})   \right\rVert_{2}  + \left\lVert  \nabla_{\bGamma} (\mathcal{L}_{n} - \bar{\mathcal{L}}_{n})(\A^{\star},\B^{\star},\bGamma^{\star})   \right\rVert_{2}  \\
&\le  \lVert \B^{\star}\rVert_{\infty} \left\rVert \mathtt{Q}_{1} \right\rVert_{2}  + n\gamma_{\max}  \left\rVert \mathtt{Q}_{2} \right\rVert_{2} +\frac{2}{2\sigma^{2}} \left\lVert  \mathtt{Q}_{4} \right\rVert_{2}.
\end{align}
Thus, combining the above bounds, we obtain
\begin{align}\label{eq:SDL_MLE_pf_bd_Q2}
S:=  \sqrt{3r} \lVert \Delta \X^{\star} \rVert_{2} + \lVert \Delta \bGamma^{\star} \rVert_{F}  \le  \sum_{i=1}^{4} c_{i} \lVert \mathtt{Q}_{i} \rVert_{2},
\end{align}
where the constants $c_{1},\dots,c_{4}>0$ are given by
\begin{align}\label{eq:c_constants_Q}
c_{1}=\left( \sqrt{3r} \lVert \B^{\star}\rVert_{\infty} + q\lVert \blambda^{\star} \rVert_{\infty} \right), \quad c_{2}=\gamma_{\max}\left( q + \sqrt{3r} \right),\quad c_{3}=q\gamma_{\max}, \quad c_{4}=\frac{2\sqrt{3}r}{2\sigma^{2}}.
\end{align}

Next, we will use concentration inequalities to argue that the right hand side in \eqref{eq:SDL_MLE_pf_bd_Q2} is small with high probability and obtain the following tail bound on $S$: %While the first three terms in \eqref{eq:SDL_MLE_pf_bd_Q2} can be made $O(1/\sqrt{n})$ with high probability, unfortunately, the last term $\mathtt{Q}_{4}$  is typically of order $O(\sqrt{p}+\sqrt{n})$ and it cannot be made small, unless the noise level $\sigma>0$ is significantly small. Nontheless, we will use concentration inequalities to obtain the following tail bound
\begin{align}\label{eq:S_tail_bound}
\P\left(S>c \sqrt{n} \log n  + 3C\sigma(\sqrt{p} + \sqrt{n}+ c\sqrt{\log n})   \right)  \le \frac{1}{n},
\end{align}
where $C>0$ is an absolute constant and $c>0$ can be written explicitly in terms of the constants we use in this proof. Recall that for a random variable $Z$, its sub-Gaussian norm, denoted as $\lVert Z \rVert_{\psi_{2}}$, is the smalleset number $K>0$ such that $\E[\exp(Z^{2}/K^{2})]\le 2$. The constant $C>$ above is the sub-gaussian norm of the standard normal variable, which can be taken as $C\le 36e/\log 2$. Using union bound with Lemmas \ref{lem:hoeffding} and \ref{lem:matrix_norm_bd}, for each $t,t>0$, we get
\begin{align}\label{eq:S_bd_pf}
&\P\left( S > (c_{1}+c_{2}+c_{3}+c_{4})t  + 3C\sigma(\sqrt{p} + \sqrt{n}+ t')   \right)  \\
&\qquad \le \left( \sum_{i=1}^{3} \P\left( \lVert \mathtt{Q}_{i} \rVert_{2}>t\right)  \right) + \P\left( \lVert n\mathtt{Q}_{4} \rVert_{2}> 3C\sigma(\sqrt{p} + \sqrt{n}+ t')    \right)  \\
&\qquad \le 2\kappa \exp\left( \frac{-t^{2} }{C_{1}^{2} \kappa^{2} n} \right) + 2p \exp\left( \frac{-t^{2}  }{ (C \sigma)^{2} p^{2} n } \right) + 2q \exp\left( \frac{-t^{2}  }{ (C \sigma')^{2}  q^{2} n } \right) + \exp(-(t')^{2}).
\end{align}
Indeed, for bounding $\P(\mathtt{Q}_{1}>t)$, we used Lemma \ref{lem:hoeffding} with sub-Gaussian norm $C_{1}=K=\gamma_{\max}/\sqrt{\log 2} $ for the bounded random vector $\dot{\h}(y_{s},\a_{s})$ (see \cite[Ex. 2.5.8]{vershynin2018high}); for $\P(\mathtt{Q}_{2}>t)$ and $\P(\mathtt{Q}_{3}>t)$, we used Lemma \ref{lem:hoeffding} with $K=C \sigma$ and $K=C \sigma'$, respectively; for the last term involving $\mathtt{Q}_{4}$, we used Lemma \ref{lem:matrix_norm_bd} with $K=C/\sigma$. Observe that in order to make the last expression in \eqref{eq:S_bd_pf} small, we will chose $t=c_{5}\sqrt{n}\log n$ and $t'=c_{5}\sqrt{\log n}$, where $c_{5}>0$ is a constant to be determined. This yields
\begin{align}
\P\left( S> c \sqrt{n} \log n  + 3C\sigma(\sqrt{p} + \sqrt{n}+ c\sqrt{\log n})   \right) \le  n^{-c_{6}},
\end{align}
where $c=c_{5}\sum_{i=1}^{4} c_{i}$ and $c_{6}>0$ is an explicit constant that grows in $c_{5}$. We assume $c_{5}>0$ is such that $c_{6}\ge 1$. This shows \eqref{eq:S_tail_bound}.

To finish, we use Theorem \ref{thm:SDL_LPGD} to deduce that with probability at least $1/n$,
\begin{align}
\lVert \bZ_{t} - \bZ^{\star}  \rVert_{F}  \le  \rho^{t}  \, \lVert  \bZ_{0} - \bZ^{\star}\rVert_{F} + \frac{\tau}{1-\rho}\left( c \sqrt{n} \log n  + 3C\sigma(\sqrt{p} + \sqrt{n}+ c\sqrt{\log n})  \right) + \frac{2\nu\tau}{1-\rho}\left( \lVert \A^{\star} \rVert_{2} + \lVert \bGamma^{\star} \rVert_{F} \right)
\end{align}
Note that $\tau<\frac{3}{2L}$ with $L = \max(2\xi,\, 2\nu+ n L^{*}) \ge n L^{*}$, so $\tau<\frac{3}{2nL^{*}}$. So this yields the desired result.

Second, suppose $\bZ^{\star}-\tau \nabla f_{\textup{SDL-filt}}(\bZ^{\star})\notin \Param$. Then we cannot direcly simplify the expression \eqref{eq:grad_mapping_compare_stationary}. In this case, we take the Frobenius norm and use non-expansiveness of the projection operator (onto convex set $\Param$):
\begin{align}\label{eq:grad_mapping_compare_stationary2}
\lVert G(\bZ^{\star}, \tau) \rVert_{F} &=  \frac{1}{\tau} \left\rVert\left[  \Pi_{\Param}\left( \bZ^{\star}-\tau\nabla \mathcal{L}_{n}(\bZ^{\star}) \right) - \Pi_{\Param}\left( \bZ^{\star}-\tau\nabla \bar{\mathcal{L}}_{n}(\bZ^{\star}) \right) \right] \right\rVert_{F}   \\
&\le \lVert \nabla \mathcal{L}_{n}(\bZ^{\star})-  \nabla \bar{\mathcal{L}}_{n}(\bZ^{\star}) \rVert_{F} \\
& \le  \lVert \Delta \X^{\star} \rVert_{F} + \lVert \Delta \bGamma^{\star} \rVert_{F}.
\end{align}
According to Remark \ref{rmk:pf_thm_LPGD}, we also have Theorem \ref{thm:CALE_LPGD} (and hence Theorem \ref{thm:SDL_LPGD}) with $\sqrt{3r} \lVert \Delta \X^{\star} \rVert_{2}$ replaced with $\lVert \Delta \X^{\star} \rVert_{F}$. Then an identical argument shows
\begin{align}\label{eq:SDL_MLE_pf_bd_Q4}
S':= \lVert \Delta \X^{\star} \rVert_{F} + \lVert \Delta \bGamma \rVert_{F} \le c_{1} \lVert \mathtt{Q}_{1} \rVert_{2}  + c_{2}\lVert \mathtt{Q}_{2} \rVert_{2} + c_{3}\lVert \mathtt{Q}_{3} \rVert_{2} + c_{4}\lVert \mathtt{Q}_{4} \rVert_{F},
\end{align}
where the constants $c_{1},\dots,c_{4}>0$ are the same as in \eqref{eq:c_constants_Q}. So we have
\begin{align}
\lVert \bZ_{t} - \bZ^{\star} \rVert_{F} \le \rho^{t}  \, \lVert  \bZ_{0} - \bZ^{\star}\rVert_{F} + \frac{\tau}{1-\rho}(S' + 2\nu (\lVert \A^{\star} \rVert_{2}+\lVert \A^{\star} \rVert_{F})).
\end{align}
Then an identical argument with the inequality $\lVert \mathtt{Q}_{4}\rVert_{F} \le \sqrt{\min(p,n)} \lVert \mathtt{Q}_{4}\rVert_{2}$ shows
\begin{align}\label{eq:S_bd_pf2}
&\P\left( S' > (c_{1}+c_{2}+c_{3}+c_{4})t  + 3C\sigma(\sqrt{p} + \sqrt{n}+ t')  \sqrt{\min(p,n)}  \right) \\
&\qquad \le \left( \sum_{i=1}^{3} \P\left( \lVert \mathtt{Q}_{i} \rVert_{2} >t\right)  \right) + \P\left( \lVert \mathtt{Q}_{4} \rVert_{2} > \frac{3C(\sqrt{p} + \sqrt{n}+ t')  }{\sigma}  \right),
\end{align}
and the assertion follows similarly as before.
\end{proof}

\begin{proof}[\textbf{Proof of Theorem \ref{thm:SDL_LPGD_STAT_feat}}]
The argument is entirely similar to the proof of Theorem \ref{thm:SDL_LPGD_STAT}. Indeed, denoting $\a_{s}=\A[:,s]+\bGamma^{T}\x_{s}'$ for $s=1,\dots,n$ and keeping the other notations the same as in the proof of Theorem \ref{thm:SDL_LPGD_STAT}, we can compute the following gradients
\begin{align}
\nabla_{\A} (\mathcal{L}_{n} - \bar{\mathcal{L}}_{n})(\A,\B,\bGamma)   &=  \begin{bmatrix}
\dot{\h}(y_{1},\a_{1}),\dots,
\dot{\h}(y_{n},\a_{n})
\end{bmatrix}  \\
\nabla_{\vect(\bGamma)} (\mathcal{L}_{n} - \bar{\mathcal{L}}_{n})(\A,\B,\bGamma) &=\left( \sum_{s=1}^{n} \dot{\h}(y_{s},\a_{s}) \otimes \x_{s}' \right) \\
\nabla_{\B}  (\mathcal{L}_{n} - \bar{\mathcal{L}}_{n})(\A,\B,\bGamma) &=\frac{2}{2\sigma^{2}} (\B^{\star}-\X_{\textup{data}}) = \frac{2}{2\sigma^{2}} [ \beps_{1},\dots,\beps_{n}] \\
\nabla_{\blambda}  (\mathcal{L}_{n} - \bar{\mathcal{L}}_{n})(\A,\B,\bGamma)  &=\frac{2}{2\sigma^{2}} \sum_{s=1}^{n}\beps_{s}' .
\end{align}
Hence repeating the same argument as before, using concentration inequalities for the following random quantities
\begin{align}
\mathtt{Q}_{1}:= \begin{bmatrix}
\dot{\h}(y_{1},\a_{1}),\dots,
\dot{\h}(y_{n},\a_{n})
\end{bmatrix}  \in \R^{p\times n} ,\quad \mathtt{Q}_{2}:=\sum_{s=1}^{n}  \beps_{s}\in \R^{p} ,\quad \mathtt{Q}_{3}:=\sum_{s=1}^{n}  \beps_{s}'\in \R^{q}, \quad \mathtt{Q}_{4}:= [ \beps_{1},\dots,\beps_{n}]  \in \R^{p\times n},
\end{align}
one can bound the size of $G(\bZ^{\star},\tau)$ with high probability. The rest of the details are omitted.
\end{proof}

\section{Proof of Theorem \ref{thm:MLE_local_consistency}}

In this section, we prove the non-asymptotic local consistency of constrained and regularized MLE, stated in Theorem \ref{thm:MLE_local_consistency}. We combine a classical approach in [21] with concentration inequalities, namely, a classical Berry-Esseen bound for deviations from standard normal distribution and  a uniform McDirmid bound (Lemma \ref{lem:uniform_McDirmid}). The former is used to control the linear term in the second-order Taylor expansion of the log-likelihood function, and the latter is used to control the second-order term. By using an $\eps$-net argument, the latter  concentration inequality can be extended to a setting where the random variables are parameterized within a compact set.

\begin{lemma}[A uniform McDirmid's inequality]
\label{lem:uniform_McDirmid}
Let $(X_{n})_{n\ge 1}$ be i.i.d. random vectors in $\R^{d}$ from a distribution $\pi$. Fix a compact parameter space $\Param\subseteq \R^{p}$ and $f_{\param}:\R^{d}\rightarrow \R$ is a bounded functional for each $\param\in \Param$ such that
\begin{align}\label{eq:McDirmid_cond}
\lVert f_{\param} - f_{\param'} \rVert_{\infty} \le L\lVert\param-\param' \rVert ,\qquad \forall \param,\param'\in \Param
\end{align}
for some constant $L>0$. Further assume that $\E_{X\sim \pi}[f_{\param}(X)]=0$ for all $\param\in \Param$. Then there exists constants $C,M>0$ such that for each $n \ge 0$, and $\eta>0$,
\begin{align}
\P\left(  \sup_{\param\in \Param} \left| \frac{1}{n}\sum_{k=1}^{n} f_{\param}(X_{k}) \right| \ge \eta \right) \le C \exp\left( p\log (1/\eta) + \frac{-\eta^{2}n}{ 2M^{2} }\right),
\end{align}
where $C'=4R(\varphi(\Omega))^{2}\tau_{\textup{min}}$.
\end{lemma}

\begin{proof}
Since $\Param\subseteq \R^{p}$ is compact, it can be covered by a finite number of $L^{2}$-balls of any given radius $\eps>0$. Let $\mathcal{U}_{\eps}$ be such an open cover using the least number of balls of radius $\eps>0$. Let $N(\eps) = |\mathcal{U}_{\eps}|$ denote the least number of such balls to cover $\Param$. Moreover, let $\textup{diam}(\Param)$ denote the diameter of $\Param$, which is finite since $\Param$ is compact. Then $\Param$ is contained in a $p$-dimensional box of side length $\textup{diam}(\Param)$. It follows that there exists a constant $K>0$, depending only on $\textup{diam}(\Param)$ and $rd$, such that
\begin{align}\label{eq:bd_N(eps)}
N(\eps)\le K\left( \frac{\textup{diam}(\Param)}{\eps} \right)^{p}.
\end{align}

Next, fix $\eta>0$, $\param\in \Param$, and $\eps>0$. Let $\param_{1},\cdots,\param_{N(\eps)}$ be the centers of balls in the open cover $\mathcal{U}_{\eps}$. Then there exists $1\le j \le N(\eps)$ such that $\lVert \param - \param_{j}\rVert <\eps$. By the hypothesis, $f_{\param}$ depends on $\param$ uniformly continuously with respect to the supremum norm. Hence there exists $\delta=\delta(\eps)>0$ such that
\begin{align}
\lVert f_{\param} - f_{\param_{j}} \rVert_{\infty} \le L\eps.
\end{align}
Denote $H_{n}(\param):= n^{-1}\sum_{k=1}^{n} f_{\param}(X_{k})$.  Then this yields, almost surely,
\begin{align}
| H_{n}(\param) - H_{n}(\param_{j}) | \le L\eps.
\end{align}
Furthermore, since $\Param$ is compact and each $f_{\param}$ is bounded, \eqref{eq:McDirmid_cond} implies that $\lVert f_{\param} \rVert_{\infty}$ is uniformly bounded in $\param$ by some constant, say, $M>0$. It follows that for each $\param\in \Param$, $H_{n}(\param)$ changes its value at most by $M$ when one of $X_{1},\dots,X_{n}$ is replaced arbitrarily. Therefore by the standard McDirmid's inequality (see, e.g., \cite[Thm. 2.9.1]{vershynin2018high}) and a union bound, with choosing $ \eps=\eta/(2L)$,
\begin{align}