\section{Introduction}

Ring lasers have been and continue to be investigated for their rich dynamics in a variety of geometries and material platforms. Numerous phenomena based on ring lasers have been numerically studied and experimentally observed, such as unidirectional and bidirectional operation, alternate oscillations, bistability and chaos~\cite{Sorel_unibistability, Sorel_opregimes, Perez_SRLdynamics, Sande_SRLdynamics, Sunada_ASEbistability, Nguimdo_chaos}. These phenomena result from the nonlinear interaction of the clockwise (CW) and counterclockwise (CCW) fields. These two fields, which in an ideal resonator do not couple to each other, give rise to the travelling-wave nature of the total field in a ring laser. This is in contrast to Fabry-Perot (FP) lasers, where the electric field exhibits a standing-wave pattern leading to spatial hole burning (SHB) of the gain medium. In lasers where the gain recovery time is fast, such as quantum cascade lasers (QCLs), SHB effects are dominant and lead to multimode instabilities~\cite{Gordon_SHB}. Although ring resonators seem to evade SHB, pure travelling waves are only possible in an ideal cavity; in a realistic cavity, even small imperfections will perturb the rotational invariance of the counter-propagating modes and pin them to a standing-wave pattern similar to FP lasers. Thus, the only realistic way to obtain travelling waves in a ring laser is by ensuring unidirectional operation where only one of the two modes, CW or CCW, is allowed to lase at a time. Unidirectional operation of ring QCLs has been shown to be necessary for the demonstration of active mode locking and mid-infrared, ring-based frequency combs~\cite{malara_extRQCL, Meng_RQCLFC}. Previous work on QCLs aimed directly at mode selection between CW and CCW modes has relied on an "S-shaped" waveguide design to favor one mode~\cite{Nshii_uniRQCL} or external cavities where the gain medium is placed within a free-space-coupled ring cavity~\cite{malara_extRQCL, Revin_activeml}.

In this work, we show robust and deterministic unidirectional operation from a monolithic racetrack QCL. Instead of a specialized waveguide design of the ring laser itself, which permanently favors one mode, we rely on an active waveguide evanescently coupled to the racetrack to perform mode selection and switching, which allows for controllable operation in \textit{both} directions. This in turn, may enable denser integration and greater functionality of mid-infrared photonic integrated circuits.

Switching between the two modes of a semiconductor ring laser has also been used to demonstrate random number generation~\cite{Sunada_ASEbistability}. The fast carrier lifetime of QCLs ($\sim$ ps) may be used in such applications to improve the speed of this operation. In addition, we show that our ring laser exhibits a region of bistability, a dynamical regime that has not yet been explored in mid-infrared QCLs. The bistable nature of the ring allows for full switching between the two modes with small input perturbations.

Our QCL operates at  $\lambda \approx$ \SI{8}{\um}, a wavelength where the fabrication of an evanescent coupler becomes challenging because the active region is buried deep in the QC layers thus requiring a complex etching process. Unlike previous work, which focuses on short wavelengths ($\lambda <$ \SI{5}{\um}), \cite{Nshii_uniRQCL, jung_integration} we fabricate an evanescent coupler that does not rely on any specialized growth or processing techniques.

Lastly, ring QCLs have recently been studied as emerging platforms for frequency combs~\cite{Meng_RQCLFC, Piccardo_phaseturb, Meng_soliton}. In these studies, light extraction is done by sensing the weak scattering off the waveguides due to bending losses. While this technique aids in minimizing backscattering between the two counter-propagating modes, it may pose a limitation on the practical use of these devices. In addition to controllable mode selection and switching, the active waveguide in our system also serves as an efficient outcoupling mechanism for the racetrack laser.

\section{Design and Fabrication}

We have fabricated a racetrack resonator coupled to a bus waveguide on a conventional $\lambda \approx$ \SI{8.8}{\um} QCL wafer based on a double-phonon resonance design \cite{Liu_qcl2phonon}. The \textit{active} bus waveguide differs from the more commonly utilized \textit{passive} bus waveguide in two key ways. First, it has intersubband transitions at the same wavelength that the racetrack emits; thus, to overcome absorption, the waveguide must be pumped to transparency. Second, due to the same transitions and depending on the pumping level, it has its own spontaneous emission which is coupled into the racetrack and is responsible for the switching properties we will discuss here.

\begin{figure}[!ht]
\centering\includegraphics[width=8.7cm]{F1schematic.pdf}
\caption{Device schematic showing the racetrack laser ridge and the evanescently-coupled waveguide with facets tilted at $17^{\circ}$. The dotted purple box shows a magnified view of the taper and s-bend employed on the ridges in order to enhance their coupling. The blue and red arrows illustrate how light is coupled in and out of the racetrack. Yellow regions represent a schematic of the top and ground metal contact pads. The bottom left panel shows the calculated mode intensity of a \SI{5}{\um} ridge waveguide. The bottom right panel offers a few combinations of gap width and coupling length leading to particular coupling strengths as calculated by CMT. The devices discussed here correspond to the last set of listed parameters.}
\label{coupler}
\end{figure}

The schematic of our system is shown in Fig.~\ref{coupler}. The two arms of the waveguide have separate metal contacts and are thus controlled independently. The coupling mechanism in this system is evanescent coupling where a small gap is deep etched between the ridges. To determine the size of the gap as well as the length of the coupler we use coupled mode theory (CMT) simulations based on the computed fundamental mode. Our CMT simulations show that significant coupling can only occur for relatively thin ridges. In our case we choose a \SI{5}{\um} ridge width because the mode is less confined laterally thereby expanding outside the waveguide core and enhancing evanescent coupling. For ease of fabrication in the following processing steps, we adiabatically taper our waveguide width to \SI{10}{\um} outside of the coupling region. This also helps to limit loss due to weak confinement and proximity of the mode to the lossy metal on the sidewalls. The detailed design of the coupler is shown in Fig.~\ref{coupler}. We have considered several combinations of gap and coupling length and their corresponding coupling strength is given in the table in Fig.~\ref{coupler}.

\begin{figure}[!ht]
\centering\includegraphics[width=8.8cm]{F2fab.pdf}
\caption{Scanning electron microscope (SEM) images of the device, focusing on the coupler. Top left: A top-view of the coupling region and taper. Top right: Cross-sectional view along the taper. The darker contour along the ridges is the insulating \ch{SiNx} layer. The coupling region is kept free of metal. Bottom left: Cross-sectional view of the ridge waveguide elsewhere in the device with a dielectric window and \ch{Ti}/\ch{Au} metal layer for electrical contact. Bottom right: Magnified view showing the etch profile and minimal, highly sub-wavelength roughness on the waveguide sidewalls.}
\label{fab}
\end{figure}

The chosen physical dimensions are compatible with conventional lithography techniques. The most challenging step in the fabrication is deep etching of the gap with an aspect ratio $\geq 5$ while still maintaining vertical sidewalls and minimal roughness. We utilize a multistep inductively coupled plasma – reactive ion etching (ICP-RIE) process based on \ch{Cl2} and \ch{BCl3}. The ICP-RIE process consists of many parameters, which we optimize for minimal sidewall roughness. Fig.~\ref{fab} shows scanning electron microscope (SEM) images of various elements of the device to especially showcase the etch characteristics. To reduce bending losses, all bend radii are designed to be larger than \SI{450}{\um}. To avoid feedback from reflection, the waveguide arms are angled $17^{\circ}$ with respect to the cleaved facet which lowers their reflectivity to $\approx 0.01$ \cite{Ahn_reflectivity, Aung_reflection}. We route the waveguides to the same side of the chip for ease of measurement, while ensuring both arms have practically the same path length to match their optical properties. Finally, an important design consideration when working with a system of two or more elements on the same wafer is electrical crosstalk due to a non-negligible ground impedance. To mitigate this issue, instead of the more conventional common, substrate-side ground used in QCLs, we use separate top ground pads for the racetrack laser and each waveguide arm.

\section{Temporal Dynamics and Mode Switching}
The racetrack laser supports two counter-propagating modes which in an ideal cavity are best described by travelling waves. Although we optimize the fabrication process, even minimal sidewall roughness will lead to backscattering and coupling between the modes. One way to preserve travelling waves in the cavity is if only one of the counter-propagating modes is allowed to lase at a time. A mechanism which allows us to select the dominant mode accurately and to efficiently switch between them is a useful tool to achieving controllable unidirectional lasing in the racetrack.
To distinguish the two modes experimentally, in our measurements we use a thermoelectrically (TE)-cooled HgCdTe (MCT) detector with a 1x\SI{1}{\mm} active area along with a pair of \ch{ZnSe} lenses of 1.5" focal length. By translating the detector across the length of the device, we are able to fully resolve the two peaks emitted from both waveguide arms, which correspond to the CW and CCW modes. All measurements are performed in pulsed mode using \SI{300}{\ns} pulses with \SI{50}{\kHz} repetition rate. The device is thermally stabilized at room temperature with a TE cooler.

We are able to detect the signal from the racetrack laser exiting the waveguide facets at ports 3 and 4 (shown in the top schematic of Fig.~\ref{switching}) even without any pump applied to the waveguide arms, although the signal is attenuated due to significant absorption along each arm. The racetrack threshold current is $\approx$ \SI{1.25}{\A} or \SI{2.89}{\kA/{\cm}^2}. Our measurements show that the racetrack laser is unidirectional; however, at each current pulse, it selects between the CW and CCW mode randomly. Under pulsed measurements, the output thus fluctuates between the two modes.

We show that a small amount of amplified spontaneous emission injected into one of the two arms favors one of the modes. In our experimental demonstration of mode selection we differentially pump the two waveguide arms taking advantage of the fact that spontaneous light emitted by each arm only couples to one of the two possible modes. Similarly, each mode will couple to a different arm of the active waveguide which allows us to distinguish them with our MCT detector. The schematic in Fig.~\ref{switching} illustrates these mode interactions with color-coded arrows. In our experiment we pump the two waveguide arms differentially by applying a variable current $I=I_{0} \pm \Delta I$ to each arm. This way the average applied current, $I_0$, at any point during our sweep remains constant. With the external MCT detector we then capture the signal at ports 3 and 4 (marked in orange in the schematic). The racetrack signal exiting the facets at 3 and 4 has interacted with the waveguide arms which provide pump-dependent absorption or amplification. An experimental measurement of this interaction is shown on the bottom right panel of Fig.~\ref{switching} which depicts the CCW output at port 3 as a function of the current applied on the left waveguide arm in the same schematic (which favors CCW operation). For up to $\approx$ 400 mA of applied current, absorption is the dominant effect; however, beyond this point, gain in the waveguide takes over and the output is amplified with increasing current. Measuring this characteristic gain curve allows us to calculate the output right as it exits the racetrack (schematically depicted by the turquoise ports 1 and 2). We use a Softplus function to fit the gain curve and then use this fit to calculate what the measured signal exiting ports 3 and 4 looks like at ports 1 and 2. This allows us to investigate and display the mode selection mechanism more clearly without the additional effects of absorption or amplification in the waveguide.

\begin{figure}[ht!]
\centering\includegraphics[width=9cm]{F3switching.png}
\caption{Top: System schematic. Bottom left: CW and CCW normalized intensity given as a function of injected current imbalance between the two waveguide arms as illustrated in the system schematic. The CW (CCW) mode shown in dotted red (blue) dominates for positive (negative) $\Delta I$. This is the signal right as it exits the racetrack, represented by the turquoise markers at ports 1 and 2 in the schematic. In solid red (blue) we show the measured probability that the mode is CW (CCW) from 15 consecutive measurements. Bottom right: CCW power measured by an MCT photodetector as it exits the waveguide (noted by an orange marker at ports 3) along with a Softplus fit.}
\label{switching}
\end{figure}

We have shown on the bottom left panel of Fig.~\ref{switching} the measured CW and CCW signals corrected using the fit discussed above as a function of the current differential between the two waveguide arms. This is the output immediately after it is outcoupled from the racetrack \textit{before} it travels through the length of the waveguide arms.
The lasing mode (CW or CCW) is only dependent on which arm has a higher pump. Following notation in the schematic of Fig.~\ref{switching}, for $\Delta I > 0$, the right arm has a higher spontaneous emission, favoring the CW mode. Conversely, for $\Delta I < 0$, the left arm's stronger emission favors the CCW mode.
For a small region around $\Delta I = 0$, the mode is chosen at random.
In the same plot we show the probability that the mode is CW (red, solid) or CCW (blue, solid) as computed from a representative set of 15 consecutive measurements for any given applied current. The intensity of the dominant mode remains relatively constant as we sweep $\Delta I$.
Any light injected into the racetrack from the waveguide arms in these measurements is due to spontaneous emission. In fact, the two waveguide arms form a FP cavity and therefore \textit{can} achieve lasing if pumped together. However, in any of these measurements, the bus waveguide is well below lasing because we only pump one arm at a time which leads to significant loss in the bus cavity. Together with the $17^{\circ}$ facets, this loss leads to a high lasing threshold ($\approx$ \SI{7.84}{\kA/{\cm}^2} or \SI{3.36}{A} if both arms are pumped simultaneously).

A mode selection mechanism based on spontaneous emission injection is consistent with our measurements. We have also considered the role that feedback of the racetrack signal from the facets back into the racetrack may play in our system; however, we have found this effect to be weak considering the combined effect of the weak reflection from the facets and a double pass through the $15\%$ coupler.

Having shown that we are able to select one of the two modes based on the current differential $\Delta I$, we now examine the effectiveness of switching between modes. We measure the amount of current imbalance needed to switch for various average currents $I_0$ applied to the waveguide arms. To this end, we assess the probability for each mode to dominate by capturing a representative train of 15 consecutive measurements from the MCT detector as shown schematically in Fig.~\ref{probability}. We only observe two possible responses from the photodetector. For example, if we are observing at the CW port, we see a high photodetector response if the mode is CW or zero response if the mode is CCW and never in between. The probability $P_{CW(CCW)}$ is then defined as the number of pulses measured in CW (CCW) mode divided by the total number of pulses, in our case 15. The results of this analysis are shown in Fig.~\ref{probability} for three different values of $I_0$. Note the crossover happens around $\Delta I \approx \SI{0}{\mA}$. It is always slightly shifted to one side, which we explain based on a slight imbalance in the lengths of the two waveguide arms. Another key detail is the span of the bistable region (where probabilities of each mode are comparable). To determine this span we fit our data using the error function and calculate the width going from $10\%$ to $90\%$ probability. For $I_0$ = \SI{620}{\mA} this region spans $\approx$ \SI{73.2}{\mA}. As $I_0$ increases to \SI{900}{\mA}, the span of this region decreases to \SI{13.8}{\mA}.

\begin{figure}[h!]
\centering\includegraphics[width=9cm]{F4probability.pdf}
\caption{Top: Schematic of a received pulse train used to measure the probability of each mode appearing. In this example, the probability $P_{CW} = P_{CCW} = 0.5$. Bottom: Probability of mode selection for the CW (red) and CCW (blue) modes. The solid line is an error function fit of the raw data (displayed as dots). The mean current increases from left to right from \SI{620}{\mA} to \SI{900}{\mA}. Required switching current, calculated as the width of the $10\%$-$90\%$ interval of each fit, decreases from \SI{73.2}{\mA} to \SI{13.8}{\mA} from left to right. The racetrack current in these measurements is held constant at \SI{1.81}{A} (\SI{4.19}{\kA/{\cm}^2}).}
\label{probability}
\end{figure}

\section{Bistability Mechanism}
Thus far we have shown that we can select the direction in which the racetrack laser emits by modulating the current differential into the waveguide arms. However, a question arises on the mechanism that enables this phenomenon. Is it necessary to continuously pump a particular arm of the waveguide if we wish to maintain the directionality of the racetrack emission or is it sustained by the racetrack dynamics itself? To answer this question we conduct the experiment shown in Fig. \ref{bistability}.
The setup is described by the upper panels of the figure. A \SI{300}{\ns} pulse of $\approx$ \SI{1.81}{\A} is applied to the racetrack. At this current, the racetrack is 1.4 times above the lasing threshold. Midway through the racetrack pulse, we apply a short current pulse of width $\approx$ \SI{30}{\ns} and current up to \SI{1.15}{A} to the waveguide arm that favors CCW operation. In the opposite arm we apply a \SI{300}{\ns} pulse of $\approx$ \SI{480}{\mA} just above the necessary current input for absorption mitigation.
Observing the photodetector response during each pulse before and after the short \SI{30}{\ns} pulse input is one way to address the question we have posed above.
On the lower left plot we show the output on the CCW side (port 3) measured on a TE-cooled MCT detector when the input CCW spike peak (in blue) is \textit{lower than} the longer CW pulse (in red). We plot 15 consecutive measurements which reveal a seemingly random selection of the lasing direction.
When the laser's CCW mode dominates, we measure a pulsed signal on the detector; conversely, when the CW mode dominates, no signal is measured, as shown in Fig.~\ref{bistability}.
On the lower-right plot we show photodetector measurements when the spike injected into the arm favoring CCW operation is \textit{higher than} the CW background injected into the other arm.
We observe that this input perturbation causes the racetrack laser to switch to the CCW mode deterministically, regardless of the initial lasing mode.
This can be seen for the section of the pulse highlighted in gray, where all of the measurements show a high signal on the detector, corresponding to the CCW mode.
Moreover, the laser remains fixed at the CCW mode, even after the perturbation that caused the switch dissipates.
This sheds some light on the mode selection mechanism.
While pumping the external waveguide arm favors one particular racetrack mode to dominate, the selected mode is self-sustaining thereafter. This demonstrates the bistable nature of the racetrack dynamics which allows for unidirectional mode selection.

\begin{figure}[h!]
\centering\includegraphics[width=9cm]{F5bistability.pdf}
\caption{Top: Input pulses into the racetrack and each waveguide arm. \textbf{(a)} The spike coupling to the CCW mode is smaller than the CW background.  \textbf{(b)} The spike coupling to the CCW mode is larger than the CW background. \textbf{(c)}, \textbf{(d)} 15 superimposed photodetector traces at port 3 corresponding respectively to pumping scenarios in (a), (b). Note the difference in the region shaded in gray between plots (c) and (d). When the short pulse exceeds pumping in the opposite arm, the CCW mode always dominates. \textbf{(e)} Device schematic along with mode interactions.
}
\label{bistability}
\end{figure}

\section{Spectral Characteristics}

\begin{figure}[!ht]
\centering\includegraphics[width=9cm]{F6spectra.pdf}
\caption{Spectral characteristics of the device. The left panel corresponds to the spectrum of the racetrack when no current is applied to the waveguide arms. For the other panels, the average applied current $I_0$ is given at the top and the spectra are given as a function of $\Delta I$ where negative (positive) $\Delta I$ favors the CW (CCW) mode.  The spectra in each panel are normalized to the maximum intensity point of that panel.
}
\label{spectra}
\end{figure}
The spectral properties of the racetrack laser and the impact that pumping of the waveguide arms has on them are summarized in Fig. \ref{spectra}. We measure the signal exiting one of the waveguide arms using a Fourier Transform Infrared Spectrometer (FTIR) with \SI{0.125}{{\cm}^{-1}} spectral resolution. The pure racetrack spectrum with no current applied to either waveguide arm is shown in the left panel and it reveals a series of longitudinal modes spanning $\approx$ \SI{14}{{\cm}^{-1}} with a measured free-spectral range (FSR) of $\approx$ \SI{0.70}{\cm^{-1}} which is an excellent match for the calculated value $1/nL$ = \SI{0.701}{\cm^{-1}}, where $L = \SI{4.318}{\mm}$ is the designed racetrack cavity length and $n = 3.3$ is the effective refractive index of the ridge waveguide. On the right panel of Fig.~\ref{spectra} we show the CW spectra as a function of current imbalance between the arms $\Delta I$ as well as average current $I_0$ which is marked at the top. In this case, a positive $\Delta I$ means higher current is applied to the waveguide arm favoring CCW operation which explains why the spectra within each panel are generally weaker for $\Delta I > 0$.

The FSR, once some current is applied to the waveguide arms, remains the same as for the pure racetrack spectrum, although there is a small red shift in the center frequency which can be explained by increased local temperature of the chip due to the applied current $2I_0$. For low values of $I_0$, despite the red shift, the spectral envelope that we measure remains similar to the racetrack spectrum without injected current. As $I_0$ increases this envelope changes significantly with only a few modes becoming dominant as can be seen for $I_0 = $ \SI{790}{\mA} and $I_0 = $ \SI{900}{\mA}. For these high values of $I_0$ we are able to recover the pure racetrack spectrum around $\Delta I = 0$ mark. This is most visible for $I_0 =$ \SI{790}{\mA} where even though only a couple of modes dominate for negative $\Delta I$, around the 0 mark, we recover the original envelope. Hence, lower values of current applied to the waveguide arms do not alter the spectral envelope of the racetrack laser; however, for higher values the envelope changes significantly. When the current differential between the arms is relatively small, the resulting spectra are also similar to the pure racetrack spectrum envelope.

\section{Conclusion}
In this paper we have shown the operation of a novel monolithic QCL-based system comprising of a racetrack laser and an evanescently-coupled bus waveguide with independently-controlled arms based on a specially-designed tapered coupler. We demonstrate deterministic and controllable selection between the two counter-propagating racetrack modes via spontaneous emission coupling from the differentially pumped waveguide arms. As little as \SI{14}{\mA} of current imbalance, which represents $\approx 1.6\%$ of the average input current into the waveguides, can trigger a full and permanenet switch between the CW and CCW modes. By further investigating the mechanism of mode selection, we have demonstrated the possibility of bistable operation of our device. Because of the bistable nature of the ring, the waveguide arms do not need constant pumping in order to maintain a particular directional mode. In fact, a small perturbation is sufficient to select the mode and the racetrack laser will maintain this mode even after the input perturbation dissipates.

Finally, we have demonstrated an efficient outcoupling mechanism by measuring up to $\approx$ \SI{4}{\mA} of optical power coupled out of the racetrack laser. These results bring us closer towards practical applications of controllable unidirectional ring or racetrack QCLs.

\paragraph{Funding} This work was supported by the Eric and Wendy Schmidt Transformative Technology Fund. S.K. acknowledges support from the Yan Huo *94 Graduate Fellowship.

\paragraph{Acknowledgments} We acknowledge the Micro and Nano Fabrication Center (MNFC) of Princeton University where fabrication was carried out.

\paragraph{Disclosures} The authors declare no conflicts of interest.

\input{arxiv.bbl}

\end{document}


\section{Introduction}\label{Intro}
\input{sections/intro}

\section{Methodology}\label{sec:methodology}
\input{sections/methodology}

\section{Efficient Parameter Estimation Algorithms}
\input{sections/alg}

\section{Asymptotic Theory}\label{sec:consistency}
\input{sections/theory}

\section{Applications}\label{sec:simulation}
\input{sections/simulation}

\section{Summary}
\input{sections/disscussion}

\section*{Acknowledgement}
The authors thank the editor, referees, and Yiyuan She for their valuable comments.
This work was partially supported by NSF grants (DMS-1924792 and CNS-1818500).

\bibliographystyle{unsrtnat}

\subsection{Related literature}

Besides regularization, another direction is to apply the low-rank constrained optimization that it has been studied for many other matrix estimation problems, including matrix sensing \citep{pmid28316458,doi:10.1137/18M1224738} and robust principal component analysis~\citep{10.5555/3291125.3309642} directly in the estimator. Compared with regularization-based methods, this direction is relatively less studied in the matrix variate regression and logistic regression settings. \citet{zhou2013tensor} applied the low-rank constraint to a tensor structure (which can be considered as a generalization of a matrix structure) and uses a block relaxation algorithm to solve the problem. However, there only exists a limited theoretical guarantee for their proposed block relaxation algorithm: consistency is only established for the unregularized estimator, and there does not exist a consistency result for the regularized estimator under the low-rank constraint. \cite{HungWang2012} proposed a maximum likelihood estimator for the model of matrix variate logistic regression based on rank-one matrix decomposition. \cite{hung2019low} applies the low-rank structure of the matrix-covariate under the same model as in this work. However, it focuses on the problem of hypothesis testing, rather than estimation as in our work.

While existing works either apply a vector-based regularization   (such as $\ell_1$ or total variation penalization \cite{li2021double}) or matrix-based  regularization such as  (such as nuclear norm penalization or low-rank constraint used in \cite{zhou2013tensor,zhou2014regularized,elsener2018robust}), the proposed estimators combine the vector-based regularization with a low-rank constraint and further improve performance. This method is consistent with the empirical observation that ``double regularization'' sometimes improves the performance, for example, the popular elastic net regularization method~\citep{zou2005regularization} uses both $\ell_1$ and $\ell_2$ regularization and achieves improved performance under certain scenarios. In addition, estimating simultaneously sparse and low-rank matrices has been studied in \cite{10.1093/imaiai/iaw012} and have applications in image processing and phase retrieval \citep{vaswani2017low}.

Our work in this paper has three main contributions. First, we propose an alternating projected gradient descent algorithm and apply it to various settings such as ordinary matrix regression, robust matrix regression, and matrix logistic regression, and the algorithm achieves the minimax estimation error when it is well-initialized. To overcome an obstacle of the nonlinear constraint in the optimization problem, we consider the constraint set in \eqref{eq:optimization} as a manifold and generalize the ``second-order Taylor expansion'' to the function defined on a manifold. Second, for the matrix variate regression problem, we investigate a method based on regularized optimization with a low-rank constraint and establish its consistency, which has not been fully studied in existing works \citep{zhou2013tensor,HungWang2012}.
In addition, we show that when the rank is small, our estimation error is smaller than the the regularization-based methods without rank constraints \citep{li2021double}, and acquires the desired statistical accuracy in a minimax sense \citep{tsybakov2008introduction,koltchinskii2011nuclear,she2021analysis}. Third, we propose a robust regularized low-rank regression model which provides accurate predictions when data contain outliers or have heavier tails than normal distributions
\citep{she2017robust}. Our consistency theorem (Theorem~\ref{prop:convergence2}) does not depend on the convexity of the loss function, and hence it applies to various choices of $l$, including popular loss function functions in robust statistics such as redescending $\psi$’s, Hampel’s loss, or Tukey’s bisquare, etc.~\citep{10.1214/aoms/1177703732,maronna2018robust,she2017robust,huang2020robust}. In addition, our theoretical analysis can be applied to general optimization problems on manifolds with bounded curvature, and thus can be considered as an important technical contribution of this paper.

\subsection{Simulation data analysis: matrix-variate regression}

In the first example, we generate a collection of 2D shapes as the true signals and the response is sampled as $y_i=\langle \bX_i,\bC^*\rangle+\langle\bz_i,\gamma^*\rangle+\epsilon_i$, where $\bX_i\in\reals^{64\times 64},\bz_i\in\reals^5$ are random with standard normal entries and $\epsilon_i$ is also standard normal. We set $\gamma^*=(1,1,1,1,1)^T$ in the following experiments. The true signal $\bC^*\in\reals^{64\times 64}$ is a 2D matrix whose entries are binary. We compare the performance of RMRE and SRRE with the various sample sizes at $n=300, 500, 700, 1000$. To tune the regularization parameter, \citet{zhou2014regularized} estimated the degree of freedom of SRRE and suggests to use Bayesian information criterion (BIC) to choose $\lambda$. In this experiment, we would rather use an independent validation set to tune the parameter $\lambda$ and a testing set to evaluate the error for both of the methods for fairness. The experiments are repeated 100 times and the mean with the standard deviation of the root mean square errors (RMSEs) of $\bC$ and $\gamma$ are present in Table~\ref{table:2d}. A comparison of the recovered signals of RMRE and SRRE at $n=500$ is shown in Fig.~\ref{fig:_triangle_circle_butterfly-robust}.

\begin{figure}[ht]
    \label{fig:2D-shapes}
    \centering
    \includegraphics[width=0.8\linewidth]{fig/2D.eps}
    \caption{Comparison of the estimators RMRE and SRRE. The first column are the true signals, the fourth column are recovered signals of RMRE without regularization. The sample size is 500. }
\end{figure}

As illustrated in Fig.~\ref{fig:_triangle_circle_butterfly-robust}, we evaluate the proposed method based on the recovered 2D shapes. First, RMRE outperforms SRRE for the square, T, and cross shapes, even without regularization ($\lambda=0$) it still estimates the signals well. RMRE and SRRE have similar performance in the triangle and circle shapes, while SRRE is slightly better for the butterfly. Second, the regularization of RMRE improved the estimation for the triangle and the circle shapes. Lastly, the recovered butterflies of both of RMRE and SRRE are barely recognized as the sample size $n=500$ is too small in this case.
The RMSEs (see  Table~\ref{table:2d}) of both $\bC$ and $\gamma$ of RMRE are much smaller than the RMSEs of SRRE with various values of $n$ for the square, T, cross shapes. As for the triangle and circle, both methods have similar performance when the sample size is small and RMRE has a better estimation as the sample size increases. RMRE performs slightly better than SRRE for the more complex shape, butterfly, when the sample size is taken as $n=1000$ and performs slightly worse when $n$ is smaller than 500. Moreover, the RMSEs of $\bC$ and $\gamma$ of RMRE are decreasing as the sample size is increasing, which verifies our theoretical analysis of the consistency in Section~\ref{sec:consistency}.

\begin{table}[ht]
\caption{Estimation errors for 2D shapes in the normal model . The means (standard deviation) of the RMSEs of $\hat \bC$ and $\hat \gamma$ with 100 repetitions of RMRE and SRRE are reported. }
\label{table:2d}
\resizebox{\linewidth}{!}{
\begin{tabular}{|c|c|c|c|c|c|c|}
\hline
Shape & Para. & Method & $n=300$ & $n=500$ & $n=700$ & $n=1000$ \\
\hline \hline
\multirow{4}{*}{Square} & \multirow{2}{*}{$\bC$}  & RMRE & 0.0144(0.0018) & 0.0105(0.0011) & 0.0093(0.0009) & 0.0086(0.0009) \\
& & SRRE & 0.1934(0.0208) & 0.0648(0.0086) & 0.0300(0.0032)   & 0.0170(0.0011) \\
\cline{3-7}
& \multirow{2}{*}{$\gamma$} & RMRE & 0.0728(0.0245) & 0.0475(0.0153) & 0.0424(0.0140)  & 0.0332(0.0103) \\
& & SRRE & 0.6465(0.2352) & 0.1686(0.0636) & 0.0815(0.0259) & 0.0429(0.0119) \\
\hline
\hline
\multirow{4}{*}{T} & $\multirow{2}{*}{$\bC$}$ & RMRE & 0.0438(0.0563) & 0.0117(0.0015) & 0.0086(0.0009) & 0.0070(0.0006)  \\
& & SRRE & 0.2068(0.0091) & 0.1388(0.0106) & 0.0762(0.0084) & 0.0328(0.0026) \\
\cline{3-7}
& $\multirow{2}{*}{$\gamma$}$ & RMRE & 0.1686(0.2022) & 0.0549(0.0172) & 0.0395(0.0129) & 0.0334(0.0107) \\
& & SRRE & 0.7381(0.2099) & 0.3768(0.1274) & 0.1844(0.0625) & 0.0716(0.0219) \\
\hline
\hline
\multirow{4}{*}{Coss} & $\multirow{2}{*}{$\bC$}$ & RMRE & 0.0368(0.0497) & 0.0113(0.0013) & 0.0083(0.0009) & 0.0065(0.0006) \\
& & SRRE & 0.1952(0.0083) & 0.1305(0.0086) & 0.0724(0.0075) & 0.0316(0.0023) \\
\cline{3-7}
& $\multirow{2}{*}{$\gamma$}$ & RMRE & 0.1386(0.1771) & 0.0547(0.0199) & 0.0398(0.0141) & 0.0321(0.0093) \\
& & SRRE & 0.7226(0.2423) & 0.3494(0.1134) & 0.1672(0.0564) & 0.0683(0.0223) \\
\hline
\end{tabular}}
\end{table}

\begin{table}[ht]
\caption{(Continued) Estimation errors for 2D shapes in the normal model. The means (standard deviation) of the RMSEs of $\hat \bC$ and $\hat \gamma$ with 100 repetitions of RMRE and SRRE are reported. }
\label{table:2d-2}
\resizebox{\linewidth}{!}{
\begin{tabular}{|c|c|c|c|c|c|c|}
\hline
Shape & Para. & Method & $n=300$ & $n=500$ & $n=700$ & $n=1000$ \\
\hline \hline
\multirow{4}{*}{Triangle} & $\multirow{2}{*}{$\bC$}$ & RMRE & 0.1777(0.0092) & 0.1036(0.0169) & 0.0615(0.0053) & 0.0546(0.0029) \\
& & SRRE & 0.1826(0.0076) & 0.1439(0.0059) & 0.1163(0.0052) & 0.0903(0.0033) \\
\cline{3-7}
& \multirow{2}{*}{$\gamma$} & RMRE & 0.5981(0.1759) & 0.2938(0.0993) & 0.1456(0.0456) & 0.1031(0.0297) \\
&  & SRRE & 0.6317(0.2171) & 0.4015(0.1238) & 0.264(0.0762)  & 0.1655(0.0587) \\
\hline \hline
\multirow{4}{*}{Circle} & \multirow{2}{*}{$\bC$} & RMRE & 0.2513(0.0101) & 0.1480(0.0195)  & 0.0529(0.0064) & 0.0413(0.0020)  \\
& & SRRE & 0.2244(0.0125) & 0.1571(0.0080)  & 0.1184(0.0058) & 0.0839(0.0039) \\
\cline{3-7}
& \multirow{2}{*}{$\gamma$} & RMRE & 0.8917(0.3156) & 0.4064(0.1590)  & 0.1303(0.0464) & 0.0880(0.0259) \\
& & SRRE & 0.8034(0.2861) & 0.4431(0.1342) & 0.2768(0.1018) & 0.1659(0.0565) \\
\hline \hline
\multirow{4}{*}{Butterfly} & \multirow{2}{*}{$\bC$} & RMRE & 0.3033(0.0036) & 0.2726(0.0067) & 0.2368(0.0106) & 0.1561(0.0124) \\
& & SRRE & 0.2894(0.0062) & 0.2602(0.0065) & 0.2312(0.0057) & 0.1954(0.0050)  \\
\cline{3-7}
& \multirow{2}{*}{$\gamma$} & RMRE & 1.0743(0.3808) & 0.7234(0.2492) & 0.5368(0.1548) & 0.3276(0.1036) \\
& & SRRE & 1.0141(0.3773) & 0.7087(0.2265) & 0.5362(0.1688) & 0.3928(0.1322) \\ \hline
\end{tabular}}
\end{table}

In the second experiment, we generate a class of synthetic signals under various ranks and sparsity levels, and apply RMRE and SRRE to these signals to compare their performance. Specifically, the matrix predictors $\bX_i\in\reals^{64\times 64}$ are matrices with each entry is standard normal and the vector predictors $\bZ_i\in\reals^5\sim\mathcal{N}(0,\mathbf{I})$ are random 5-dimensional vectors. The number of sample is chosen as $n=500$. Moreover, we let $\gamma^*=(1,1,1,1,1)^T$ and $\bC^*=\bC_1\bC_2^T$, where $\bC_d\in\{0,1\}^{64\times r},d=1,2$. Each entry of $\bC_d$ is taken as 1 with probability $\sqrt{1-(1-s)^{1/r}}$ and 0 otherwise. The parameter $r$ controls the rank of $\bC$ and $s$ controls its sparsity. In the experiment, we choose $r=1,5,10,20$ and $s=0.01,0.05,0.1,0.2,0.5$. We consider a  normal model where the response $y_i$ is generated as $\langle \bX_i,\bC^*\rangle+\gamma^{*T}\bZ_i+\epsilon_i$, where $\epsilon_i$ is standard normal.

We evaluate the performance of RMRE and SRRE with respect to prediction. We trained both methods on a train set with $500$ samples and evaluated the prediction errors on an independent test set with sample size $500$. The prediction error is defined as the misclassification rate of the responses in the logistic model and is defined as RMSE of the responses in the normal model. The regularization parameter $\lambda$ is chosen as the prediction error on the test set is minimized and this approach is applied for both methods for fairness. We reported the mean and the standard deviations of the prediction errors with 100 repetitions in Table~\ref{table:matrix-syn-3} for the normal model. The results for RMSEs of the matrix coefficient $\bC^*$ are reported in Table~\ref{table:matrix-syn-4}.
Compared with SRRE, our proposed method has good estimation results for the rank-1 matrices in both normal and logistic models. As the sparsity decreases, our method outperforms SRRE for matrix of various ranks. When the sparsity parameter $s$ is less or equal than $0.2$, our method at least has a similar performance both in RMSEs in estimating $\bC^*$ and prediction in the responses. The only exception case is for large sparsity matrices ($s=0.5$ and $r\geq 5$), SRRE has smaller estimation errors.

\begin{table}[ht]
\caption{Estimation result for synthetic data in the logistic model. Reported are the mean and standard deviation (in the parenthesis) of the rate of prediction error of $\hat{Y}_i$ out of 100 data replications. RMRE is the proposed method, regularized matrix-variate regression estimator, and SRRE is the spectral regularized regression estimator. }
\label{table:matrix-syn-1}
\begin{center}
\resizebox{\linewidth}{!}{
\begin{tabular}{|cc|c|c|c|c|}
\hline
Sparsity & Method & \multicolumn{4}{c|}{Rank} \\ \cline{3-6}
& & $r=1$ & $r=5$ & $r=10$ & $r=20$  \\ \hline \hline
$1\%$    & RMRE   & 0.1178(0.0218) & 0.2504(0.0449) & 0.2999(0.0359) & 0.3346(0.0346) \\
         & SRRE   & 0.2351(0.0298) & 0.3454(0.0302) & 0.3641(0.0252) & 0.3755(0.0260)  \\ \hline \hline
$5\%$    & RMRE   & 0.0994(0.0178) & 0.3460(0.0338) & 0.3899(0.0255) & 0.3994(0.0247) \\
         & SRRE   & 0.2034(0.0226) & 0.3688(0.0287) & 0.3952(0.0230) & 0.4053(0.0256) \\ \hline \hline
$10\%$   & RMRE   & 0.1192(0.0234) & 0.3744(0.0317) & 0.4059(0.0233) & 0.4098(0.0230) \\
         & SRRE   & 0.2012(0.0260) & 0.3704(0.0258) & 0.3986(0.0227) & 0.4092(0.0226) \\ \hline \hline
$20\%$   & RMRE   & 0.1347(0.0216) & 0.3799(0.0259) & 0.4060(0.0279) & 0.4146(0.0204) \\
         & SRRE   & 0.1978(0.0257) & 0.3474(0.0303) & 0.3828(0.0297) & 0.4008(0.0249)\\ \hline \hline
$50\%$   & RMRE   & 0.1293(0.0192) & 0.3605(0.0262) & 0.3935(0.0226) & 0.4107(0.0206) \\
         & SRRE   & 0.1915(0.0253) & 0.2844(0.0316) & 0.3145(0.0351) & 0.3467(0.0306)  \\ \hline
\end{tabular}}
\end{center}
\end{table}

\begin{table}[ht]
\caption{Estimation result for synthetic data in the logistic model. Reported are the mean and standard deviation (in the parenthesis) of the RMSEs of $\hat{\bC}$ out of 100 data replications. RMRE is the proposed method, regularized matrix-variate regression estimator, and SRRE is the spectral regularized regression estimator. }
\label{table:matrix-syn-2}
\begin{center}
\resizebox{\linewidth}{!}{
\begin{tabular}{|cc|c|c|c|c|}
\hline
Sparsity & Method & \multicolumn{4}{c|}{Rank} \\ \cline{3-6}
 & & $r=1$ & $r=5$ & $r=10$ & $r=20$ \\ \hline \hline
$1\%$    & RMRE   & 0.0697(0.0250) & 0.0842(0.0199) & 0.0921(0.0188) & 0.0954(0.0140) \\
         & SRRE   & 0.0860(0.0256) & 0.0915(0.0189) & 0.0970(0.0183) & 0.0984(0.0137) \\ \hline \hline
$5\%$    & RMRE   & 0.1920(0.0328) & 0.2225(0.0286) & 0.2238(0.0264) & 0.2235(0.0223) \\
         & SRRE   & 0.2090(0.0317) & 0.2251(0.0280) & 0.2248(0.0264) & 0.2238(0.0222)\\ \hline \hline
$10\%$   & RMRE   & 0.2817(0.0431) & 0.3301(0.0367) & 0.3315(0.0273) & 0.3345(0.0301) \\
         & SRRE   & 0.3019(0.0464) & 0.3327(0.0376) & 0.3321(0.0272) & 0.3350(0.0301) \\ \hline \hline
$20\%$   & RMRE   & 0.4090(0.0447) & 0.4972(0.0411) & 0.5159(0.0398) & 0.5105(0.0411) \\
         & SRRE   & 0.4421(0.0466) & 0.5021(0.0420) & 0.5176(0.0400) & 0.5110(0.0411)\\ \hline \hline
$50\%$   & RMRE   & 0.6568(0.0384) & 0.9765(0.0672) & 1.0326(0.0727) & 1.0458(0.0635)\\
         & SRRE   & 0.6932(0.0388) & 0.9855(0.0678) & 1.0343(0.0725) & 1.0464(0.0634)\\ \hline
\end{tabular}}
\end{center}
\end{table}

\begin{table}[ht]
\caption{Estimation result for synthetic data in the normal model. Reported are the mean and standard deviation (in the parenthesis) of the RMSEs of responses $\hat{Y}_i$ out of 100 data replications. RMRE is the proposed method,  regularized matrix-variate regression estimator, and SRRE is the spectral regularized regression estimator. }
\label{table:matrix-syn-3}
\begin{center}
\resizebox{\linewidth}{!}{
\begin{tabular}{|cc|c|c|c|c|}
\hline
Sparsity & Method & \multicolumn{4}{c|}{Rank}\\ \cline{3-6}
 & & $r=1$ & $r=5$ & $r=10$ & $r=20$ \\ \hline \hline
$1\%$    & RMRE   & 1.0664(0.0433) & 1.5320(0.5631) & 2.5858(1.0295) & 5.0129(1.1507) \\
         & SRRE   & 2.1023(0.2856) & 5.1372(1.1028) & 5.6477(1.0383) & 5.9977(0.9132) \\ \hline \hline
$5\%$    & RMRE   & 1.0925(0.0483) & 4.4520(3.5651) & 11.7671(1.9385) & 13.2728(1.5353)\\
         & SRRE   & 3.1843(0.4954) & 12.0552(1.4539) & 12.8586(1.4045) & 13.5944(1.3929)\\ \hline \hline
$10\%$   & RMRE   & 1.1874(0.0907) & 11.8509(5.5520) & 19.1758(2.0377) & 20.0153(1.6428)\\
         & SRRE   & 4.2707(0.9060) & 17.4163(1.8125) & 19.1120(1.8045) & 19.7686(1.5268)\\ \hline \hline
$20\%$   & RMRE   & 1.5175(0.1847) & 25.2633(3.1545) & 29.7988(2.5752) & 31.1058(2.6314)\\
         & SRRE   & 5.7842(1.2161) & 25.2321(2.1572) & 27.9915(2.2142) & 29.8019(2.4316)\\ \hline \hline
$50\%$   & RMRE   & 2.9755(0.3597) & 52.1691(4.8101) & 59.9110(4.4410) & 63.3346(4.5564)\\
         & SRRE   & 8.3832(1.6463) & 41.1190(2.9915) & 49.2962(3.3167) & 54.4517(3.3581)\\ \hline
\end{tabular}}
\end{center}
\end{table}

\begin{table}[ht]
\caption{Estimation result for synthetic data in the normal model. Reported are the mean and standard deviation (in the parenthesis) of the RMSEs of $\hat{\bC}$ out of 100 data replications. RMRE is the proposed method, regularized matrix variate regression estimator, and SRRE is the spectral regularized regression estimator. }
\label{table:matrix-syn-4}
\begin{center}
\resizebox{\linewidth}{!}{
\begin{tabular}{|cc|c|c|c|c|}
\hline
Sparsity & Method & \multicolumn{4}{c|}{Rank} \\ \cline{3-6}
&  & $r=1$ & $r=5$ & $r=10$ & $r=20$ \\ \hline \hline
$1\%$    & RMRE   & 0.0053(0.0009) & 0.0174(0.0102) & 0.0368(0.0170) & 0.0765(0.0179)\\
         & SRRE   & 0.0285(0.0047) & 0.0784(0.0174) & 0.0860(0.0160) & 0.0923(0.0135)\\ \hline \hline
$5\%$    & RMRE   & 0.0067(0.0010) & 0.0661(0.0569) & 0.1815(0.0303) & 0.2049(0.0222)\\
         & SRRE   & 0.0471(0.0079) & 0.1878(0.0225) & 0.1989(0.0215) & 0.2101(0.0196)\\ \hline \hline
$10\%$   & RMRE   & 0.0099(0.0024) & 0.1836(0.0876) & 0.2975(0.0316) & 0.3114(0.0265)\\
         & SRRE   & 0.0648(0.0139) & 0.2708(0.0275) & 0.2969(0.0262) & 0.3080(0.0244)\\ \hline \hline
$20\%$   & RMRE   & 0.0177(0.0035) & 0.3955(0.0475) & 0.4640(0.0370) & 0.4833(0.0369)\\
         & SRRE   & 0.0888(0.0189) & 0.3932(0.0294) & 0.4354(0.0317) & 0.4627(0.0331)\\ \hline \hline
$50\%$   & RMRE   & 0.0433(0.0055) & 0.8153(0.0672) & 0.9296(0.0593) & 0.9862(0.0630)\\
         & SRRE   & 0.1284(0.0244) & 0.6426(0.0440) & 0.7645(0.0436) & 0.8468(0.0459)\\ \hline
\end{tabular}}
\end{center}
\end{table}

\subsection{Simulation data: matrix-variate logistic regression}

We consider a  logistic regression model where the response $y_i$ is $1$ with probability $\frac{1}{1+e^{-\langle \bX_i,\bC^*\rangle-\gamma^{*T}\bZ_i}}$ and 0 otherwise.
We evaluate the performance of our method (RMRE) and SRRE, and report the mean and the standard deviations of the predication errors under 100 replications in Table~\ref{table:matrix-syn-1} for the logistic model. The results for RMSEs of the matrix coefficient $\bC$ are reported in Table~\ref{table:matrix-syn-2}.

\subsection{Simulation data analysis: robust matrix-variate regression}

We compare the proposed robust matrix variate regression with the robust low-rank matrix estimator (RLRME)~\citep{elsener2018robust}, which can be considered a robust version of SRRE by replacing the quadratic loss with the Huber loss.
The value of the tuning parameter is suggested as $\lambda=2\sqrt{\log(m+q)/q/n}$ in \citet{elsener2018robust}. For a fair comparison, we instead use the same method to tune the parameter in the following simulations. We trained both of Robust RMRE and RLRME on a train set with certain number of samples and evaluated the prediction errors on an independent validation set with the sample number of samples.  The results are in Tables 1 to 5. The responses in both of the train set and the validation set are generated by
$y=\langle \bX,\bC^*\rangle+\gamma^{*T}\bZ+(1-p)\epsilon+p\epsilon'$, where $\epsilon'\sim\mathcal{N}(0,10^4)$ is the contamination with probability $p$. The results are in Tables 6 and 7. The value of $\gamma^*$ is set to $(1,1,1,1,1)^T$ for all simulations in this section.

We use RMSEs of $\bC$ and $\gamma$ with the various sample sizes to compare the proposed method, robust regularized matrix-variate regression estimator (Robust RMRE), with RLRME, which requires a certain number of samples to make an accurate estimation. Yet, our proposed method does not have such a requirement and has smaller errors even the sample size is small. In this experiment, the true signal $\bC^*\in\reals^{30\times 30}$ is a square shape and is defined as $\bC^*_{ij}=1$ if $10\leq i,j\leq 20$ and $\bC^*_{ij}=0$ otherwise.
The noise is set as $\epsilon\sim\mathcal{N}(0,1)$ and the contamination probability is $p=0.05$ and $p=0.1$.
Fig.~\ref{fig:msre-sample-size} shows the RMSEs of the estimated $\bC$ and $\gamma$ against the sample size with the corrupted responses. The response has the corruption with probability 0.05 in the left figure of Fig.~\ref{fig:msre-sample-size} and this probability is taken as 0.1 in the right figure. In both figures, RLRME has large errors when the sample size is smaller than 1000 and our method (RobustRMRE) always has small errors even the sample size is also small. This phenomenon can be explained by Theorem~\ref{thm:asymptotics}, where it is shown that Robust RMRE has a better consistency rate than RLRME when $r$ is smaller than $m$ and $q$.

We further compare the robust estimators RobustRMRE and RLRME with non-robust estimator RMRE for various 2D shapes (Triangle, Circle, Butterfly). The same 2D shape signals in \citet{zhou2013tensor} and \citet{zhou2014regularized} are used here but the matrix size is changed to $40\times 40$. The sample size is taken as 1000 and the responses are corrupted with probability $0.05$. Fig.~\ref{fig:_triangle_circle_butterfly-robust} exhibits that RobustRMRE can recover better signals than RLRME while the non-robust estimator has a poor estimation when the outliers are present. The numerical results are reported in Table~\ref{table:robust_gaussian_noise} with the various sample sizes. The responses of the training and validation sets are corrupted with probability $0.05$. We conducted 100 repetitions of the simulation and reported the mean and the standard deviation of the RMSEs. Overall, the errors of our method, RobustRMRE, are much smaller than RLRME and the non-robust estimator, RMRE. RLRME is robust to the outliers for most of the shapes except for some simple shapes such as the square. Both RobustRMRE and RLRME have smaller errors as the sample size is increasing. An additional comparison is made when the noise distribution is chosen as the standard Cauchy distribution. The RMSEs of $\mathbf{C}$ and $\gamma$ against the sample size with Cauchy noise are shown in Fig.~\ref{fig:msre-sample-size-cauchy}. The experimental results reported in Table~\ref{table:robust_cauchy} also indicate that our method RobustRMRE outperforms than RLRME when the noise distribution is heavy-tailed.

\begin{figure}[ht]
   \centering
   \includegraphics[width=0.45\linewidth]{fig/p0.05.eps}
   \includegraphics[width=0.45\linewidth]{fig/p0.1.eps}
   \caption{Comparison of the proposed robust estimator (Robust RMRE) with RLRME as the number of samples increases. The outliers are generated with probability 0.05 (left) and 0.1 (right). Here $C_1:=\hat{\bC}^{RobustRMRE}, C_2:=\hat{\bC}^{RLRME}, \theta_1:=\hat{\gamma}^{RobustRMRE},\theta_2:=\hat{\gamma}^{RLRME}$ denote the estimated values of $\bC$ and $\gamma$. The curves are obtained out of 10 replications. }
   \label{fig:msre-sample-size}
\end{figure}

\begin{figure}[ht]
   \centering
   \includegraphics[width=0.6\linewidth]{fig/cauchy-noise.eps}
   \caption{Comparison of the proposed robust estimator (Robust RMRE) with RLRME as the number of samples increases. The noise follows the standard Cauchy distribution. Here $C_1:=\hat{\bC}^{RobustRMRE}, C_2:=\hat{\bC}^{RLRME}, \theta_1:=\hat{\gamma}^{RobustRMRE},\theta_2:=\hat{\gamma}^{RLRME}$ denote the estimated values of $\bC$ and $\gamma$. The curves are obtained out of 10 replications. }
   \label{fig:msre-sample-size-cauchy}
\end{figure}

\begin{figure}[ht]
  \centering
  \includegraphics[width=\linewidth]{fig/_triangle_circle_butterfly-robust.eps}
  \caption{Comparison of the robust estimators RobustRMRE and RLRME with non-robust estimator. The sample size is 1000. }
  \label{fig:_triangle_circle_butterfly-robust}
\end{figure}

\begin{table}[ht]
\caption{Estimation errors for 2D shapes in the normal model with corruption probability 0.1. The means (standard deviation) of the RMSEs of $\hat \bC$ and $\hat \gamma$ with 100 repetitions are reported. }
\label{table:robust_gaussian_noise}
\resizebox{\linewidth}{!}{
\begin{tabular}{|c|c|c|c|c|c|c|}
\hline
Shape & Para. & Method & $n=300$ & $n=500$ & $n=700$ & $n=1000$\\
\hline \hline
\multirow{6}{*}{Triangle} & \multirow{3}{*}{$\bC$} & RobustRMRE & 0.1440(0.0150) & 0.0885(0.0168) & 0.0600(0.0077) & 0.0622(0.0062) \\
& & RLRME & 0.2172(0.0116) & 0.1990(0.0092)  & 0.1727(0.0064) & 0.1427(0.0057) \\
& & RMRE & 0.6123(0.1349) & 0.8563(0.1645) & 0.7337(0.1248) & 0.5503(0.0650) \\
\cline{2-7}
& \multirow{3}{*}{$\gamma$} & RobustRMRE & 0.3419(0.1251) & 0.1851(0.0711) & 0.1106(0.0419) & 0.0899(0.0347) \\
& &  RLRME  & 0.6011(0.1842) & 0.4174(0.1220)  & 0.3609(0.1090)  & 0.2900(0.0810) \\
& & RMRE & 1.8787(0.7706) & 1.7311(0.7102) & 1.3096(0.4453) & 0.9755(0.3208) \\
\hline \hline
\multirow{6}{*}{Circle} &\multirow{3}{*}{$\bC$}    & RobustRMRE & 0.1924(0.0159) & 0.1075(0.0185) & 0.0606(0.0098) & 0.0403(0.0032) \\
& & RLRME & 0.2763(0.0082) & 0.2490(0.0099)  & 0.2166(0.0083) & 0.1796(0.0084) \\
& & RMRE & 0.6470(0.1362)  & 0.9113(0.1528) & 0.7231(0.1106) & 0.5448(0.0644) \\
\cline{2-7}
& \multirow{3}{*}{$\gamma$} & RobustRMRE & 0.4352(0.1286) & 0.2192(0.0761) & 0.1079(0.0416) & 0.0670(0.0227)  \\
& & RLRME      & 0.6440(0.1726)  & 0.4870(0.1541)  & 0.3972(0.1222) & 0.3192(0.0920)  \\
& & RMRE       & 1.9324(0.7693) & 1.8282(0.7066) & 1.2984(0.4442) & 0.9762(0.3591)  \\
\hline \hline
\multirow{6}{*}{Butterfly} &\multirow{3}{*}{$\bC$} & RobustRMRE & 0.2766(0.0095) & 0.2309(0.0103) & 0.1632(0.0209) & 0.1149(0.0137) \\
&  & RLRME & 0.2967(0.0094) & 0.2726(0.0112) & 0.2410(0.0081)  & 0.1996(0.0076) \\
& & RMRE & 0.4828(0.0621) & 0.6332(0.0928) & 0.8299(0.1162) & 0.9296(0.1105) \\
\cline{2-7}
& \multirow{3}{*}{$\gamma$} & RobustRMRE & 0.5529(0.1547) & 0.4366(0.1175) & 0.2683(0.0966) & 0.1761(0.0574) \\
& & RLRME & 0.6431(0.1860)  & 0.5522(0.1642) & 0.4051(0.1167) & 0.3338(0.1027) \\
& & RMRE & 1.5754(0.5818) & 1.4130(0.5340)   & 1.3980(0.5066)  & 1.3328(0.4371) \\ \hline
\end{tabular}}
\end{table}

\begin{table}[ht]
\caption{Estimation errors for 2D shapes in the normal model as the noise follows standard Cauchy distribution. The means (standard deviation) of the RMSEs of $\hat \bC$ and $\hat \gamma$ with 100 repetitions are reported. }
\label{table:robust_cauchy}
\resizebox{\linewidth}{!}{
\begin{tabular}{|c|c|c|c|c|c|c|}
\hline
Shape & Para. & Method & $n=300$ & $n=500$ & $n=700$ & $n=1000$ \\
\hline \hline
\multirow{6}{*}{Triangle} & \multirow{3}{*}{C} & RobustRMRE & 0.1643(0.0128) & 0.1171(0.0130)  & 0.0685(0.0143)  & 0.0607(0.0085) \\
& & RMRE & 1.7292(4.3080) & 3.1972(6.7697) & 4.6344(24.463) & 3.0650(9.2047)\\
& & RLRME & 0.2219(0.0052)  & 0.2011(0.0054) & 0.1837(0.0057) & 0.1580(0.0048)\\
\cline{2-7}
& \multirow{3}{*}{$\gamma$} & RobustRMRE & 0.3992(0.1189) & 0.2530(0.0990) & 0.1486(0.0623) & 0.1182(0.0456) \\
& & RMRE & 5.5389(12.511) & 6.4430(12.667) & 7.9979(42.893) & 5.0391(15.485) \\
& & RLRME & 0.5862(0.1694) & 0.4667(0.1285) & 0.4025(0.1099) & 0.3076(0.1011) \\
\hline \hline
\multirow{6}{*}{Circle} & \multirow{3}{*}{C} & RobustRMRE & 0.2158(0.0162) & 0.1985(0.0154) & 0.0923(0.0142) & 0.0611(0.0257) \\
& & RMRE & 1.4626(5.3714) & 5.0158(14.751) & 2.0511(2.8052) & 12.075(98.787) \\
& & RLRME & 0.2847(0.0052) & 0.2593(0.0072) & 0.2320(0.0078) & 0.1987(0.0118) \\
\cline{2-7}
& \multirow{3}{*}{$\gamma$} & RobustRMRE & 0.5134(0.1821) & 0.3874(0.1319) & 0.1917(0.0660)  & 0.1263(0.0908) \\
& & RMRE & 3.7961(9.1388) & 9.0904(21.233) & 3.9441(5.6502) & 15.248(115.30) \\
& & RLRME & 0.6609(0.1794) & 0.5335(0.1282) & 0.4126(0.1229) & 0.3333(0.1068) \\
\hline \hline
\multirow{6}{*}{Butterfly} & \multirow{3}{*}{C} & RobustRMRE & 0.2834(0.0081) & 0.2950(0.0255) & 0.2035(0.0132) & 0.1264(0.0128) \\
& & RMRE & 1.1195(2.8638) & 2.2435(5.7978) & 9.5673(44.680) & 3.1074(4.5313) \\
& & RLRME & 0.2946(0.0047) & 0.2706(0.0058) & 0.2461(0.0070) & 0.2143(0.0061) \\
\cline{2-7}
& \multirow{3}{*}{$\gamma$} & RobustRMRE & 0.5635(0.1931) & 0.5175(0.1635) & 0.3371(0.1005) & 0.1942(0.0755) \\
& & RMRE & 4.3562(13.464) & 4.7767(10.976) & 13.589(49.213) & 4.7824(10.478)\\
& & RLRME & 0.6330(0.1829) & 0.4982(0.1299) & 0.4534(0.1238) & 0.3440(0.1035)\\
\hline
\end{tabular}}
\end{table}

\subsection{Real Data Analysis}

We compared RMRE, SRRE, and SDNCMV on the electroencephalography (EEG) data of alcoholism \citep{zhang1995event}.
There are 77 individuals with alcoholism and 44 individuals as the control (non-alcoholic). The subjects are exposed to a stimulus and the voltage values were measured from 64 channels of electrodes placed on the subject's scalp for 256 time points and 120 trials. We average these sampling data out of 120 trials and collect 122 matrices of size $256\times 64$ of the individuals. The responses of these individuals are taken as either 0 (non-alcoholic) or 1 (alcoholic).
The classical linear model deals with vector-valued predictors, which will lead to poor performance if we simply vectorize the matrix sampling data into a vector. For instance, the dimension of the sampling unit is $256\times64=16,384$ but the sample size is only 122. Moreover, vectorization destroys the wealth of structural information inherently possessed in the matrix data~\citep{zhou2014regularized}.

We apply RMRE, SRRE, and SDNCMV to the data and evaluate the prediction error on the testing set via the cross-validation. We used the same modified cross-validation method in \cite{zhou2014regularized} to achieve a fairer way in tuning the regularization parameter. Specifically, we divided the full data into a training and a testing
sample using $k$-fold cross-validation.
Then for the training data, we further employed a 5-fold cross-validation
to tune the shrinkage parameter $\lambda$. We then applied the tuned model that is fully based
on the training data now to the testing data and evaluated the misclassification error rate
for testing. We report the average performance on the random partition of the training and testing sets in cross validation with $100$ repetitions.

\begin{table}[t!]
\caption{Misclassification rates (standard deviation in parentheses)  of RMRE, SRRE, and SDNCMV for the EEG data}
\label{table:eeg-data}
\resizebox{\linewidth}{!}{\begin{tabular}{|c|c|c|c|c|}
\hline
Method & Leave-one-out CV & 5-fold CV & 10-fold  CV & 20-fold CV \\ \hline \hline
\textbf{RMRE} & 0.2131 & 0.2298 (0.0255) & 0.2255 (0.0214) & 0.2186 (0.0104) \\
SRRE & 0.2131 & 0.2332 (0.0225) & 0.2228 (0.0188) & 0.2171 (0.0122) \\
SDNCMV &  0.2459 & 0.2348 (0.0064) & 0.2349 (0.0071) & 0.2610 (0.0086) \\
\hline
\end{tabular}}
\end{table}

We also apply RMRE, SRRE, and SDNCMV to classify the IEEE leucorrhoea microscopic images \citep{hao2019automatic},
which consist of 6 categories: Erythrocytes (Ery, also known as red blood cells), Leukocytes (Leu, also known as white blood cells), Molds, Epithelial Cells (Epi), Pyocytes (Pyo) and Trichomonads. Examples of these images are present in Figure~\ref{fig:leu-sample} where Trichomonads are not found in this dataset. To evaluate the  performance, we sample $120$ images from the categories of Leu and Pyo ($60$ images each) and downsample the resolution of each image to $32$ by $32$.
The comparisons on the EEG data and the leucorrhoea data for leave-one-out, 5, 10, and 20-fold cross-validation are reported in Tables~\ref{table:eeg-data} and \ref{table:leucorrhoea-data}, where the mean misclassification rates and its standard deviation over 100 runs are recorded (the standard deviation of leave-one-out CV is zero as the training/testing partition and the algorithms are deterministic).

\begin{figure}[ht]
  \centering
  \includegraphics[width=0.19\linewidth,height=0.12\textheight]{fig/1-Ery.eps}
  \includegraphics[width=0.19\linewidth,height=0.12\textheight]{fig/2-Leu.eps}
  \includegraphics[width=0.19\linewidth,height=0.12\textheight]{fig/3-Epi.eps}
  \includegraphics[width=0.19\linewidth,height=0.12\textheight]{fig/4-Mid.eps}
  \includegraphics[width=0.19\linewidth,height=0.12\textheight]{fig/5-Pyo.eps}
  \caption{Microscopic images for each category. The types of images are shown as Erythrocytes (Ery, \nth{1}), Leukocytes (Leu, \nth{2}), Epithelial Cells (Epi, \nth{3}), Mildew (Mid, \nth{4}), Pyocytes (Pyo, \nth{5})}
  \label{fig:leu-sample}
\end{figure}

\begin{table}[t!]
\caption{Misclassification rates (standard deviation in parentheses)  of RMRE, SRRE, and SDNCMV for the leucorrhoea data}
\label{table:leucorrhoea-data}
\resizebox{\linewidth}{!}{
\begin{tabular}{|c|c|c|c|c|}
\hline
Method & Leave-one-out CV & 5-fold CV & 10-fold  CV & 20-fold CV \\ \hline \hline
\textbf{RMRE} & 0.2750 & 0.2293 (0.0213) & 0.2537 (0.0204) & 0.2528 (0.0172)\\
SRRE & 0.2583 & 0.2605 (0.0251) & 0.2593 (0.0149) & 0.2582 (0.0161) \\
SDNCMV & 0.2583 & 0.2583 (0.0081) & 0.2461 (0.0106) & 0.2653 (0.0076) \\
\hline
\end{tabular}}
\end{table}
\subsection{Consistency of the proposed estimator}\label{sec:mainresults}
We first present a few conditions for establishing the model consistency of the proposed estimator in \eqref{eq:optimization}.

\textbf{Assumption 1}: There exists a positive constant $C_1>0$ such that
\[
\frac{1}{n}\Big\|\sum_{i=1}^n \mathrm{vec}(\bX_i,\bz_i)\mathrm{vec}(\bX_i,\bz_i)^T\Big\|\leq C_1,
\]
where $\mathrm{vec}(\bX_i,\bz_i)\in\reals^{qm+p}$ is a vector consisting of the elements in $\bX_i\in\reals^{m\times q}$ and $\bz_i\in\reals^p$.

\textbf{Assumption 2}: For $\bH(\bC,\gamma)\in\reals^{(qm+p)\times (qm+p)}$  defined by
\begin{equation}\label{eq:hessian}
\bH(\bC,\gamma)=\sum_{i=1}^nw_{2,i}\mathrm{vec}(\bX_i,\bz_i)\mathrm{vec}(\bX_i,\bz_i)^T,
\end{equation}
where
\begin{align*}
w_{2,i}=\begin{cases}2,\,\,\text{for the ordinary matrix variate regression model,}\\
I\Big(|y_i-\langle\bX_i,\bC\rangle-\bz_i^T\gamma|<\alpha\Big),\,\,\text{for the robust matrix variate regression model,}\\
e^{\langle\bX_i,\bC\rangle+\bz_i^T\gamma}/(1+e^{\langle\bX_i,\bC\rangle+\bz_i^T\gamma})^2,\,\,\text{for the logistic matrix variate regression model},\end{cases}
\end{align*}
and $\frac{1}{n}\bH(\bC,\gamma)$ is positive definite with eigenvalues bounded from below for all  $(\bC,\gamma)$ in a neighborhood of $(\bC^*,\gamma^*)$. Specifically, there exists positive constants $C_2>0$ and $c_0\leq \sigma_r(\bC^*)/2$ such that
$
\frac{1}{n}\lambda_{\min}\Big(\bH(\bC,\gamma)\Big)\geq C_
$ and for all $(\bC,\gamma)$ such that $$\dist((\bC,\gamma),(\bC^*,\gamma^*))=\sqrt{\|\bC-\bC^*\|_F^2+\|\gamma-\gamma^*\|^2}\leq c_0.$$

Assumptions 1 and 2 can be considered as the generalized version of Restricted Isometry condition in \cite{doi:10.1137/070697835} and are comparable to  \cite[Conditions 1,2,5,6]{li2021double}. In particular, Assumption 1 ensures that the sensing vectors $\mathrm{vec}(\bX_i,\bz_i)\in\reals^{qm+p}$ are not the same or in a similar direction, and it is satisfied if the sensing vectors are sampled from a distribution that is relatively uniformly distributed among all directions. Assumption 2 ensures that the Hessian matrix of $\sum_{i=1}^n l(y_i,\langle \Xb_i,\Cb\rangle+\gamma^T\zb_i)$, $\bH$, is not a singular matrix.

\begin{remark}
Assumptions 1 and 2 are reasonable when $n>O(pm+q)$ and $(\bX_i,\bz_i)$ is sampled from a reasonable distribution that does not concentrate around certain directions. For example, if $\mathrm{vec}(\bX_i,\bz_i)$ are i.i.d. sampled from a distribution of $N(0,\bI)$, then the standard concentration of measure results~\citep{wainwright_2019} imply that for the matrix variate regression model, $\bH=2\sum_{i=1}^n \mathrm{vec}(\bX_i,\bz_i)\mathrm{vec}(\bX_i,\bz_i)^T$, and the standard concentration of measure results~\citep{wainwright_2019} imply that Assumption 1 holds with a high probability when $n=O(qm+p)$ and Assumption 2 holds for the regression model as well since $\bH=2\sum_{i=1}^n \mathrm{vec}(\bX_i,\bz_i)\mathrm{vec}(\bX_i,\bz_i)^T$. Assumption 2 holds for the logistic regression model as long as  $a_i=\langle\bX_i,\bC\rangle+\bz_i^T\gamma$ is bounded above for most indices $i$ as $w_{2,i}\geq 1/(2e^{a_i})$, and holds for the robust regression model if  $\alpha$ is not too small such that there are sufficient samples (e.g., more than $2(qm+p)$ samples) such that $|y_i-\langle\bX_i,\bC\rangle-\bz_i^T\gamma|<\alpha$.
\end{remark}

\textbf{Assumption 3} [Assumption on the noise for the matrix variate regression model and the robust matrix variate regression model]: For the matrix variate regression model, the error $\epsilon_i$'s in Equation \eqref{eq:model1} follow an independent and identically distributed (i.i.d.) zero-mean and sub-Gaussian distributions with zero mean and variance one, i.e., $\Expect (\epsilon_i )= 0$ and $\Var (\epsilon_i )=\sigma^2$ to ensure that the distribution of outliers is limited as this model is not designed to handle outliers. For the robust matrix variate regression model, the errors $\epsilon_i$'s follow i.i.d. zero-mean and possibly heavy-tailed distributions.

With Assumptions 1-3, we establish the consistency of the proposed  estimator \eqref{eq:optimization} in Theorem~\ref{thm:asymptotics} that shows for large $n$, our estimator converges to the true underlying solution $(\bC^*,\gamma^*)$. This result holds for general penalty $P$ and regularization parameter $\lambda$.  A sketch of the proof is presented in Section~\ref{sec:skecth} and the detailed proof is deferred to the Appendix

\begin{theorem}[Consistency of the proposed estimator]\label{thm:asymptotics}
Under Assumptions 1-3, let
\begin{align*}
\sigma_0=\begin{cases}\sigma,\,\,\text{for the matrix variate regression model}\\1,\,\,\text{for the robust matrix variate regression model}\\1,\,\,\text{for the logistic matrix variate regression model},\end{cases}
\end{align*}
then for
\begin{align}
n\geq C \left(6\frac{C_1\sigma_0\sqrt{(r(qm+p)}}{C_2}+3\sqrt{\frac{\lambda P(\bC^*,\gamma^*)}{C_2}}\right)^2\Bigg/c_0^
 \label{eq:assumption_asymptotics}
\end{align}
and for all $t\geq 2$, we have the following upper bound on the estimation error of the proposed estimator \eqref{eq:optimization} with probability at least $1-C\exp(-Cn)-C\exp(-Ct(r(q+m)+p))$:
\begin{equation}\label{eq:bound_consistency}
\dist((\hat{\bC},\hat{\gamma}),(\bC^*,\gamma^*))
\leq 6\frac{C_1 t \sigma_0\sqrt{(r(q+m)+p)} }{C_2\sqrt{n}}+\sqrt{\frac{6\lambda P(\bC^*,\gamma^*)}{C_2n}}.
\end{equation}
\end{theorem}
\begin{comment}
{\color{red}
\begin{definition}[convergence rate??]
  For any arbitrary $\varepsilon>0$,
\begin{equation*}
Pr\left(\sqrt{\|\hat{\bC}-\bC^*\|_F^2+\|\hat{\gamma}-\gamma^*\|^2} > \varepsilon r_n
\right)\rightarrow 0\mbox{~~as~~} n\rightarrow \infty.
\end{equation*}
\end{definition}
\noindent Note that the consistency is at a rate $r_n$ such that $r_n \rightarrow 0$ as $n \rightarrow \infty$ implies consistency.
}
\end{comment}
\subsection{Discussion of Theorem~\ref{thm:asymptotics}}
\hfill

\textbf{Comparison of the convergence rate with existing works.}
The main result, Inequality \eqref{eq:bound_consistency}, shows that the estimation error is in the order of \begin{equation}\label{eq:error_order}O_P\left(\sqrt{\frac{r(q+m)+p}{n}}+\sqrt{\frac{\lambda P(\bC^*,\gamma^*)}{n}}\right).\end{equation} In comparison, the rate in the existing work \citep{li2021double} is $$O_P\left(\sqrt{\frac{qm+p}{n}}+\sqrt{\frac{\lambda P(\bC^*,\gamma^*)}{n}}\right).$$
Our result improves the factor $mq + p$ to $r(q+m)+p$, which is a significant improvement when the rank $r$ is smaller than $q$ and $m$. This improvement can be explained by the fact that by fixing the rank constraint, our estimator has $r(q+m)+p$ parameters which are fewer than the $qm+p$ parameters in the estimator of \cite{li2021double}. With fewer parameters to estimate, our estimator achieves a better convergence rate, and this improvement is also clear from our simulation experiments. $O_P$ and $o_P$ represent the standard big and small O notation in probability.

\textbf{Comparison with minimax rate.}  We derive our minimax analysis in Theorem \ref{prop:convergence2}, which shows that under the setting that $\bC^*$ is low-rank, the minimax rate of estimating $(\bC,\gamma)$ is in the order of $O_P\left(\sqrt{\frac{r(q+m)+p}{n}}\right)$. In comparison, our rate in Theorem~\ref{thm:asymptotics}  achieves this rate. In literature, \citet[Remark 10]{luo2020recursive} shows that the minimax rate is in the order of $O_P\left(\sqrt{\frac{r(q+m)}{n}}\right)$ when the vector covariates $\gamma$ does not exist and $p=0$ and can be considered as a special case of our result. In the future, we plan to derive an improved minimax rate of our model under the setting that $\bC^*$ is both sparse and low-rank based on the technique of \cite{she2021analysis}, and show that our estimator with a well-chosen $\lambda\neq 0$ and $P(\bC,\gamma)=\|\bC\|_1$ achieves  the improved rate.

\textbf{MLE argument when $\lambda=0$.} In addition, the adaption of the usual arguments of MLE to manifold optimization can be applied when $\lambda=0$, in which we may show that $\sqrt{n}(\tilde{\bC}-\bC,\tilde{\gamma}-\gamma)$ converges to a Gaussian-like distribution (the result is similar to \citet[Theorem 3.1]{HungWang2012}), where the covariance/Fisher information matrix can be obtained following the technique of \cite{le1990maximum} and \cite{6418045}. However, we skip the rigorous statement here as the we use $\lambda>0$ in practice.

\textbf{A generic choice of $P$.}  In this work, we leave the choice of the penalty $P$ open for the theoretical analysis. However, there are natural choices of $P$ in certain applications. For example, when we have the prior knowledge that $\bC^*$ is sparse, then we may let $P(\bC,\gamma)=\|\bC\|_1$.

\textbf{Choice of rank $r$} The choice of $r$ has a large impact on the solution. In fact, the objective value of \eqref{eq:optimization} is nonincreasing as $r$ increases, since the set of matrix of rank $r_2$ contains the set of matrices of rank $r_1$ when $r_2>r_1$. We expect that a choice of $r$ can be obtained using either cross validation, or the elbow method \citep{choi2017selecting} that chooses the rank such that a larger rank doesn't give much smaller objective value.

\subsection{Convergence of the proposed algorithm}
In this section, we show that the output of Algorithm~\ref{alg:matrix-reg-log} is linearly convergent and sufficiently close to the ground truth with a good initialization and an additional assumption on the second order information of the ``Hessian matrix'' $\bH$ in Assumption 2.

\textbf{Assumption 4}: There exists a positive constant $C_3>0$ such that $\frac{1}{n}\Big\|\bH(\bC,\gamma)\Big\|\leq C_3$ for all $(\bC,\gamma)$ such that $\dist((\bC,\gamma),(\bC^*,\gamma^*))=\sqrt{\|\bC-\bC^*\|_F^2+\|\gamma-\gamma^*\|^2}\leq c_0.$
Combining it with Assumption 2, it suggests that $\bH$ is a well-conditioned matrix. Hence, Assumption 4 can be considered as the generalized version of the Restricted Isometry condition in \cite{doi:10.1137/070697835} and comparable to  \citet[Conditions 2,6]{li2021double}.

Assumptions 4 is reasonable when $n>O(pm+q)$ and $(\bX_i,\bz_i)$ is sampled from a reasonable distribution that does not concentrate around certain directions. For example,  if $\mathrm{vec}(\bX_i,\bz_i)$ are i.i.d. sampled from a distribution of $N(0,\bI)$, then the standard concentration of measure results~\citep{wainwright_2019} imply that Assumption 2 holds for all three models with a high probability, since $w_{2,i}$ are bounded above.
With Assumption 4, we have the following result showing the convergence of Algorithm~\ref{alg:matrix-reg-log}, with its proof deferred to the Appendix. It shows that with a good initialization, all accumulation points have estimation errors converging to zero as $n\rightarrow\infty$.

\begin{theorem}[Convergence of the proposed algorithm]\label{prop:convergence2}
Under Assumptions 1-4, let $c_0'=\min\left(c_0,\frac{C_2}{4C_TC_3}\right)$,
the initialization $(\bC^{(0)},\gamma^{(0)})$ is ``good'' in the sense that it is close to $(\bC^*,\gamma^*)$, the objective value $F(\bC^{(0)},\gamma^{(0)})$ is not large:
\begin{align}\nonumber
\dist(({\bC}^{(0)},{\gamma}^{(0)}),(\bC^*,\gamma^*))&\leq c_0',\\
F(\bC^{(0)},\gamma^{(0)})-\sum_{i=1}^n l(y_i,\langle \Xb_i,\Cb^*\rangle+\gamma^{*T}\zb_i)&\leq \frac{C_2n}{2}c_0'^2,\label{eq:convergence_init}
\end{align}
and $n$ is large: $\sqrt{n}>\frac{4C_TC_1\sigma_0\sqrt{qm+p} }{C_2}.$
Then for all $t\geq 2$, with a probability at least $1-C\exp(-Cn)-C\exp(-t(r(q+m)+p))$, all accumulation points of the iterates $\{(\bC^{(\iter)}, \gamma^{(\iter)})\}_{\iter\geq 1}$, denoted by $(\tilde{\bC},\tilde{\gamma})$, have small estimation errors:
\begin{equation}\label{eq:stationary}
\dist((\tilde{\bC},\tilde{\gamma}),(\bC^*,\gamma^*))\leq \frac{4C_1t\sigma_0\sqrt{r(q+m)+p} }{C_2\sqrt{n}}+\frac{4\lambda C_{partial}}{nC_2},
\end{equation}
 where $C_{partial}$ is the upper bound of the magnitude of the derivative of the penalty function within a neighborhood of the ground truth:
$$C_{partial}=\max_{\{\dist((\bC,\gamma),(\bC^*,\gamma^*))\leq c_0'\}}\sqrt{\Big\|\frac{\partial}{\partial \bC}P(\bC,\gamma)\Big\|_F^2+\Big\|\frac{\partial}{\partial \gamma}P(\bC,\gamma)\Big\|^2}.$$
\end{theorem}
\begin{remark}
Unlike Theorem~\ref{thm:asymptotics}, Theorem~\ref{prop:convergence2} does not depend on the convexity of the loss function $l$. As a result, it applies to various choices of $l$, including popular loss function functions in robust statistics such as redescending $\psi$’s, Hampel’s loss, or Tukey’s bisquare, etc.~\citep{10.1214/aoms/1177703732,maronna2018robust,she2017robust,huang2020robust} that can detect outliers with moderate or high leverages.
\end{remark}

\subsection{Minimax rate}\label{sec::minimax}
In this section, we prove the minimax lower bounds
showing that the rates attained by our estimator are optimal, based on the arguments in \cite{tsybakov2008introduction}. We consider the
class of parameters
\[
\mathcal{A}(r,a)=\{(\bC,\gamma)\in\reals^{q\times p}\times \reals^p: \mathrm{rank}(\bC)\leq r, \|\bC\|_F^2+\|\gamma\|^2\leq 1\}.
\]
We will need the following assumptions:\\
\textbf{Assumption 5:} There exists $C_{upper}>0$ such that for all $(\bC,\gamma)\in \mathcal{A}(r,a)$,
\[
 \sum_{i=1}^n\Big(\langle\bX_i,\bC\rangle+\gamma^T\bz_i\Big)^2\leq  C_{upper}n\|\mathrm{vec}(\bC,\gamma)\|^2.
\]
\textbf{Assumption 6:} There exists $c_{\epsilon}>0$ such that for all $x\in\reals$,
$
KL(P_{\epsilon,0},P_{\epsilon,x})\leq c_{\epsilon}x^2,
$
where $P_{\epsilon,x}$ represent the distribution of  $\epsilon_i+x$ for the regression model and the robust regression model, and $P_{\epsilon,x}$ represent the Bernoulli distribution with parameter $x$ for the logistic regression model. In addition, $KL$ represents the Kullback-Leibler divergence: For distributions $P$ and $Q$ of a continuous random variable with probability density functions $p(x)$ and $q(x)$, it is defined to be the integral $KL(P,Q)=\int_{x}p(x)\log(p(x)/q(x))\di x$.

Assumption 5 is less restrictive of Assumption 1 as it only need to be true for all $(\bC,\gamma)\in\mathcal{A}(r,a)$, and Assumption 6 holds under the logistic regression model, and under the regression and robust regression models, Assumption 6 holds for  zero-mean, symmetric distributions with tails decaying not faster than Gaussian, including Gaussian distribution, exponential distribution, Cauchy distribution, and Student's t distribution.
Now we are ready to state our main result on the minimax rate. Combining it with convergence rate in Theorem~\ref{thm:asymptotics}, it implies that the rate of our estimator is optimal.
\begin{theorem}\label{thm:minimax}
Assuming that $n\geq r(m+q)+p$, then for any $0<\beta<\frac{2^{(r(m+q)+p)/4}}{1+2^{(r(m+q)+p)/4}}$, there exists $c_0>0$ depending on $\beta, c_{\epsilon}$ and $C_{upper}$ such that
\[
\inf_{\hat{\bC},\hat{\gamma}}\sup_{(\bC^*,\gamma^*)\in \mathcal{A}(r,a)}\Pr\left(\dist\Big((\bC^*,\gamma^*),(\hat{\bC},\hat{\gamma})\Big)\geq c_0 \sqrt{\frac{r(m+q)+p}{n}} \right)\geq \beta.
\]
\end{theorem}
\begin{proof}
WLOG assume that $m\geq q$. Let $\theta=(\bC,\gamma)$, and
\[
\mathcal{C}=\Bigg\{(\bC,\gamma):\bC=[\bC',0]\,\, \text{where}\,\, \bC'\in\reals^{m\times r}\,\, \text{and}\,\, 0\in\reals^{m\times (q-r)},  \bC'_{ij}\in\{s,-s\}, \gamma_k\in\{s,-s\}\Bigg\},
\]
where $s=c_0\sqrt{\frac{r(m+q)+p}{n}}$. Then we have $|\mathcal{C}|=2^{rm+p}$, and for any $(\bC_1,\gamma_1), (\bC_2,\gamma_2)\in\mathcal{C}$, $\dist((\bC_1,\gamma_1), (\bC_2,\gamma_2))\geq 2s$. In addition, for any $(\bC,\gamma)\in\mathcal{C}$ and $P_0$ represents the model when $\bC=0$ and $\gamma=0$,
\[
K(P_{0},P_{(\bC,\gamma)})\leq c_{\epsilon}\sum_{i=1}^n\Big(\langle\bX_i,\bC\rangle+\gamma^T\bz_i\Big)^2\leq c_{\epsilon}C_{upper}(\|\bC\|_F^2+\|\gamma\|^2)\leq c_{\epsilon}C_{upper} (r(m+q)+p)s^2.
\]
Applying \cite[Theorem 2.5]{tsybakov2008introduction} and note that $\log|\mathcal{C}|=2^{rm+p}=(rm+p)\log 2$, we may choose $
\alpha={2c_{\epsilon}C_{upper}}{s^2}/\log 2$, and $\alpha$ can be sufficiently small by choosing $c_0$ to be small. The rest of the proof following applying \cite[Theorem 2.5]{tsybakov2008introduction}.
\end{proof}
\subsection{Sketch of the Proof of Theorem~\ref{thm:asymptotics}}\label{sec:skecth}
We start with the intuition of the proof with a function $f: \reals^p\rightarrow\reals$. To show that $f$ has a local minimizer around $\bx^*$, it is sufficient to show that the gradient $\nabla f(\bx^*)\approx 0$ and the Hessian matrix of $f(\bx)$, $\bH(\bx)$, is positive definite with eigenvalues strictly larger than some constant $c>0$. The intuition of the proof follows from the Taylor expansion that \begin{align}\label{eq:taylor}f(\bx)\approx& f(\bx^*)+(\bx-\bx^*)^T\nabla f(\bx^*) + \frac{1}{2}(\bx-\bx^*)^T H(\bx^*)(\bx-\bx^*)\\\geq& f(\bx^*)+(\bx-\bx^*)^T\nabla f(\bx^*) + \frac{c}{2}\|\bx-\bx^*\|^2.\nonumber\end{align} As a result, there is local minimizer in the neighbor of $\bx^*$ with radius $\frac{\|\nabla f(\bx^*)\|}{c}$, i.e., $B\left(\bx^*,\frac{\|\nabla f(\bx^*)\|}{c}\right)$.
To extend this proof to \eqref{eq:optimization}, the main obstacle is the nonlinear constraint in the optimization problem. To address this issue, we consider the constraint set in \eqref{eq:optimization} as a manifold and generalize the ``second-order Taylor expansion'' in \eqref{eq:taylor} to the function defined on a manifold. With this generalized Taylor expansion, a similar strategy can be applied to prove that the minimizer of \eqref{eq:optimization} is close to $(\bC^*,\gamma^*)$.

To analyze functions defined on manifolds, we introduce a few additional notations. We assume a manifold $\calM\subset\reals^p$ and a function $f:\reals^p\rightarrow \reals$, and investigate $f(\bx)$ for $\bx\in B(\bx^*,r)\cap \calM$, i.e., a local neighborhood of $\bx^*$ on the manifold $\calM$. We denote the first and second derivatives of $f(\bx)$ by $\nabla f(\bx)\in\reals^{p}$ and $\bH(\bx)\in\reals^{p\times p}$ respectively, the tangent plane of $\calM$ at $\bx^*$ by $T_{\bx^*}(\calM)$, and let $\Pi_{T_{\bx^*}(\calM)}$ and $\Pi_{T_{\bx^*}(\calM),\perp}$ be the projectors to $T_{\bx^*}(\calM)$ and its orthogonal subspace respectively. These definitions are visualized in Figure~\ref{fig:manifold}.

\begin{figure}[ht]
  \centering
  \includegraphics[width=0.7\linewidth]{fig/manifold.eps}
  \caption{A visualization of the manifold $\calM$, two points $\bx,\bx^*\in\calM$, the tangent space $T_{\bx^*}(\calM)$, and the projectors $\Pi_{T_{\bx^*}(\calM)}$ and $\Pi_{T_{\bx^*}(\calM),\perp}$. }
  \label{fig:manifold}
\end{figure}

Then, we say that a manifold $\calM$ is curved with parameter $(c_0,C_T)$ at $\bx^*\in\calM$, if for any $\bx\in B(\bx^*,c_0)\cap\calM$, we have $\|\Pi_{T_{\bx^*}(\calM),\perp}(\bx-\bx^*)\|\leq C_T\|\Pi_{T_{\bx^*}(\calM)}(\bx-\bx^*)\|^2$. Intuitively, it means that the projection of $\bx-\bx^*$ to the tangent space $T_{\bx^*}(\calM)$ has a larger magnitude than the projection to the orthogonal subspace of the tangent space (see Figure~\ref{fig:manifold}).We remark that larger $C_T$ means that the manifold $\calM$ is more ``curved'' around $\bx^*$. Then,  Lemma~\ref{lemma:general_manifold} establishes the lower bound of $f$ based on the local properties such as the first and second derivatives of $f$ at $\bx^*$, the tangent space of $\calM$ around $\bx^*$, and the curvature parameters $(c_0,C_T)$.
\begin{lemma}\label{lemma:general_manifold}
Consider a $d$-dimensional manifold $\calM\subset\reals^p$ and a function $f: \reals^P\rightarrow\reals$, define
$
C_{H,1}= \min_{\bx\in B(\bx^*,c_0)} \lambda_{\min}(\bH(\bx)), $
and assume that $\calM$ is curved with parameter $(c_0,C_T)$ at $\bx^*$ and $4C_{H,1}\geq C_T\|\Pi_{T_{\bx^*}(\calM),\perp}\nabla f(\bx^*)\|$, then we have the following lower bound for any $\bx\in B(\bx^*,c_0)\cap\calM$:
\[
f(\bx)-f(\bx^*)\geq \frac{1}{2}b^2C_{H,1}- b \|\Pi_{T_{\bx^*}(\calM)}\nabla f(\bx^*)\|-C_Tb^2\|\Pi_{T_{\bx^*}(\calM),\perp}\nabla f(\bx^*)\|,
\]
where $b=\|\bx-\bx^*\|$.
\end{lemma}

Lemma~\ref{lemma:general_manifold} can be viewed as a generalization of Inequality \eqref{eq:taylor}: when $\calM=\reals^p$, then we have $C_T=0$, $\Pi_{T,\perp}=\emptyset$ and as a result, $\|\Pi_{T_{\bx^*}(\calM),\perp}\nabla f(\bx^*)\|=0$.
To apply Lemma~\ref{lemma:general_manifold} to our problem, we need to estimate the parameters  $c_0, C_T, C_{H,1},$  $\|\Pi_{T_{\bx^*}(\calM)}\nabla f(\bx^*)\|$, and $\|\Pi_{T_{\bx^*}(\calM),\perp}\nabla f(\bx^*)\|$ in the statement of Lemma~\ref{lemma:general_manifold}. In particular, $(c_0,C_T)$ depends on the manifold $\calM$ used in the optimization problem \eqref{eq:optimization}, which is \begin{equation}\calM=\{(\bC,\gamma)\in\reals^{q\times m}\times \reals^p: \rank(\bC)=r\}\label{eq:manifold}.\end{equation} By treating $\calM$ as the product space of $\reals^p$ and the manifold of low-rank matrices $\{\bC\in\reals^{q\times m}: \rank(\bC)=r\}$ and following the  tangent space  of the set of low-rank matrices in the literature \citep{Absil2015,10.5555/3291125.3309642}, we obtain the following lemma of ``curvedness'' parameters $(c_0,C_T)$ of $\calM$ at $(\bC^*,\gamma^*)$.

\begin{lemma}\label{lemma:curvature}
The manifold $\calM$ defined in \eqref{eq:manifold} is curved with parameter $(c_0,C_T)$ at $(\bC^*,\gamma^*)$, for any $c_0\leq \sigma_{r}(\bC^*)/2$ and $C_T=2/\sigma_{r}(\bC^*)$, where $\sigma_{r}(\bC^*)$ represents the $r$-th (i.e., the smallest) singular value of $\bC^*$.
\end{lemma}
In addition, the parameter $C_{H,1}$ in Lemma~\ref{lemma:general_manifold} can be estimated from Assumption 2, and the derivatives $\|\Pi_{T_{\bx^*}(\calM)}\nabla f(\bx^*)\|$ and $\|\Pi_{T_{\bx^*}(\calM),\perp}\nabla f(\bx^*)\|$ in Lemma~\ref{lemma:general_manifold} can be estimated from Assumptions 1 and 3. Then the proof of Theorem~\ref{thm:asymptotics} follows from  Lemma~\ref{lemma:general_manifold} and the intuition introduced at the beginning of this section, with technical details deferred to the Appendix.

\subsection{Discussion of Theorem~\ref{prop:convergence2}}\label{sec:convergence2}
\hfill

\textbf{Convergence.} The key result in this section, \eqref{eq:stationary}, shows that the algorithm achieves a similar estimation error as the result in the consistency result in Inequality \eqref{eq:bound_consistency}, with the term $\sqrt{\frac{6\lambda P(\bC^*,\gamma^*)}{C_2n}}$ replaced with $\frac{4\lambda C_{partial}}{nC_2}$. For many common penalty functions, $C_{partial}$ is bounded. For example, for the $\ell_1$ loss function that $P(\bC,\gamma)=\|\bC\|_1$ used in our simulations, we have $C_{partial}\leq \sqrt{qm}$.

\textbf{Condition on initialization in Inequality \eqref{eq:convergence_init}.} Theorem~\ref{prop:convergence2} makes two assumptions on the initialization in Inequality \eqref{eq:convergence_init}, where the LHS of the first condition is the distance between $(\bC^{(0)},\gamma^{(0)})$ and the true solution $(\bC^*,\gamma^*)$. While the second condition is more complicated, but it is satisfies when $\dist((\bC^{(0)},\gamma^{(0)}),(\bC^*,\gamma^*))=o_P(1)$ and $\lambda=o(1)$ as $n\rightarrow\infty$: note that the LHS is
\begin{align*}
&F(\bC^{(0)},\gamma^{(0)})-\sum_{i=1}^n l(y_i,\langle \Xb_i,\Cb^*\rangle+\gamma^{*T}\zb_i)\\=&\sum_{i=1}^n l(y_i,\langle \Xb_i,\Cb^{(0)}\rangle+\gamma^{(0)T}\zb_i)-\sum_{i=1}^n l(y_i,\langle \Xb_i,\Cb^*\rangle+\gamma^{*T}\zb_i)+\lambda P(\bC^{(0)},\gamma^{(0)}),
\end{align*}
and Assumption 4 implies that
\[
\sum_{i=1}^n l(y_i,\langle \Xb_i,\Cb^{(0)}\rangle+\gamma^{(0)T}\zb_i)-\sum_{i=1}^n l(y_i,\langle \Xb_i,\Cb^*\rangle\leq nC_3 \dist((\bC^{(0)},\gamma^{(0)}),(\bC^*,\gamma^*))^2.
\]

\textbf{Initialization.} Theorem~\ref{prop:convergence2} requires an initialization that is within a neighborhood of the true solution of radius $O(1)$. While empirically we initialize $\bC^{(0)}$ and $\gamma^{(0)}$ as a zero matrix and a zero vector, it is possible to obtain initialization with theoretical guarantees. For example, for the regression setting, we may let the initialization to be the solution of the standard regression problem, i.e., the solution of $\argmin_{\bC,\gamma}\sum_{i=1}^n (y_i-\langle \Xb_i,\Cb\rangle-\gamma^T\zb_i)^2$, which has an initial estimation error in the order of $O_P\left(\sqrt{\frac{qm+p}{n}}\right)$. For the robust regression and logistic regression settings, we may solve $\argmin_{\bC,\gamma}\sum_{i=1}^n l(y_i,\langle \Xb_i,\Cb\rangle+\gamma^T\zb_i)$ as well. These are convex problems, and the standard asymptotic analysis shows that the estimation errors would converge to zero as $n\rightarrow\infty$, that is, the conditions in the inequality \eqref{eq:convergence_init} will be satisfied as $n\rightarrow\infty$.

\textbf{Convergence rate.} Unfortunately, it is difficult to obtain a general result of the convergence rate to $(\hat{\bC},\hat{\gamma})$ without assumptions on the penalty $P$. However, our proof implies that Algorithm~\ref{alg:matrix-reg-log} converges linearly to a neighborhood around the optimal solution $(\bC^*,\gamma^*)$, and empirically Algorithm~\ref{alg:matrix-reg-log} converges quickly and the convergence rate are linear in our simulations.