text:
\section{Applications to GLMs -- \emph{tossing-coin-sketch} method} \label{sec:GLM}

Consider the problem of training a generalized linear model
\begin{equation} \label{eq:P}
w^* = \arg\min_{w \in \R^d} P(w) \, \eqdef \, \frac{1}{n}\sum_{i=1}^n \phi_i(a_i^\top w) + \frac{\lambda}{2}\|w\|^2,
\end{equation}
where $\phi_i:\R \rightarrow \R^+$ is a convex and continuously twice differentiable loss function, $a_i \in \R^d$ are data samples and $w \in \R^d$ is the parameter to optimize. As the objective function is strongly convex, the unique minimizer satisfies $\nabla P(w) = 0$, that is
\begin{align}\label{eq:nablaP}
\nabla P(w) &= \frac{1}{n}\sum_{i=1}^n\phi_i'(a_i^\top w)a_i + \lambda w = 0.
\end{align}
Let
$
\Phi(w) \;\eqdef\;
\begin{bmatrix}
\phi_1'(a_1^\top w) \ \cdots \ \phi_n'(a_n^\top w)
\end{bmatrix}^\top \;\in \;\R^{n}
$
and
$
\mA  \;\eqdef\;
\begin{bmatrix}
a_1 \ \cdots \ a_n
\end{bmatrix} \;\in \; \R^{d \times n}.
$
By introducing auxiliary variables $\alpha_i \in \R$ s.t.\  $\alpha_i \eqdef -\phi_i'(a_i^\top w)$, we can re-write~\eqref{eq:nablaP} as
\begin{equation}
w \; = \; \frac{1}{\lambda n}\mA\alpha,
\quad \quad \mbox{and } \quad \quad
\alpha \; = \; -\Phi(w). \label{eq:fixedPhi}
\end{equation}
Note $x = [\alpha; w] \in \R^{n+d}$.
The objective of finding the minimum of \eqref{eq:P} is now equivalent to finding \emph{zeros} for the function

\begin{equation}\label{eq:F}
F(x) \; = \; F(\alpha; w) \; \eqdef \;
\begin{bmatrix}
\frac{1}{\lambda n}\mA\alpha - w \\
\alpha + \Phi(w)
\end{bmatrix},
\end{equation}
where $F:\R^{n+d} \rightarrow \R^{n+d}$. \ruiline{Our objective now becomes solving $F(x) = 0$ with $p = m = n+d$.} For this, we will use a variant of the \SNR. The advantage in representing~\eqref{eq:nablaP} as the nonlinear system~\eqref{eq:F} is that we now have one row per data point (see the second equation in~\eqref{eq:fixedPhi}). This allows us to use sketching to \emph{subsample} the data.

Since the function $F$ has a block structure, we will use a structured sketching matrix which we refer to as a \emph{Tossing-coin-sketch}. But first, we need the following definition of a block sketch.

\begin{definition}[($n,\tau$)--block sketch]
Let $B_n \subset \{1,\ldots, n\}$ be a subset of size $\tau$ uniformly sampling at random. We say that $\mS \in \R^{n\times \tau}$ is a $(n,\tau)$--block sketch if $\mS = \mI_{B_n}$ where $\mI_{B_n}$ denotes the column concatenation of the columns of the identity matrix $\mI_n \in \R^{n\times n}$ whose indices are in  $B_n$.
\end{definition}

Our Tossing-coin-sketch is a sketch that alternates between two blocks depending on the result of a coin toss.

\begin{definition}[Tossing-coin-sketch]\label{def:TCS}
Let $\mS_d \in \R^{d \times \tau_d}$ and $\mS_n \in \R^{n \times \tau_n}$ be a $(d,\tau_d)$--block sketch and a  $(n,\tau_n)$--block sketch, respectively. Let $b \in (0; 1)$.

Now each time we sample $\mS$, we ``toss a coin'' to determine the structure of $\mS  \in \R^{(d+n)\times(\tau_d+\tau_n)}$. That is, $\mS = \begin{bmatrix}
\mS_d & 0 \\
0 & 0
\end{bmatrix}$ with probability $1-b$ and $\mS = \begin{bmatrix}
0 & 0 \\
0 & \mS_n
\end{bmatrix}$ with probability $b$.
\end{definition}

By applying the \SNR method with a tossing-coin-sketch for solving~\eqref{eq:F}, we arrive at an efficient method for solving~\eqref{eq:P} that we call the \emph{TCS} method.

By using a tossing-coin-sketch, we can alternate between solving a linear system based on the first $d$ rows of~\eqref{eq:F} and a nonlinear system based on the last $n$ rows of~\eqref{eq:F}.

\ruilinex{TCS is inspired by the first-order stochastic dual ascent methods~\cite{SDCA,Shalev-Shwartz2015,Qu2015b}. Indeed, eq.~\eqref{eq:fixedPhi} can be seen as primal-dual systems with primal variables $w$ and dual variables $\alpha$. Stochastic dual ascent methods are efficient to solve~\eqref{eq:fixedPhi}. At each iteration, they update alternatively the primal and the dual variables $w$ and $\alpha$ with the first-order informations. Thus, by sketching alternatively the primal and the dual systems and updating accordingly with the Newton-type steps, TCS's updates can be seen as the second-order stochastic dual ascent methods.}

We show in the next section that the TCS method verifies~\eqref{eq:DFTinEmHass}.

Using sketch sizes s.t.\ $\tau_n \ll n$,

the TCS method has the same cost as \SGD in the case $d \ll n$.
\ruiline{The low computational cost per iteration is thus another advantage of the TCS method.}
See App.~\ref{sec:detailexp} the cost per iteration analysis.
For a detailed derivation of the TCS method, see App.~\ref{sec:equivalencealgorithms} and a detailed implementation in Algorithm~\ref{algo:tau-TCS} in the appendix.

\subsection{The condition~\texorpdfstring{\eqref{eq:DFTinEmHass}}{} in the case of TCS method}
\label{sec:kerTCS}

In this section, we show that the TCS method verifies~\eqref{eq:DFTinEmHass} through Lemma~\ref{lem:reformulationnew} in the following.
\ruiline{
\begin{lemma} \label{lem:adaptTCS}
Consider the function $F$ defined in \eqref{eq:F} and the tossing-coin-sketch $\mS$ defined in Definition~\ref{def:TCS}, then~\eqref{eq:DFTinEmHass} holds.
\end{lemma}}

\begin{proof}
First, we show that $\E{\mS\mS^\top}$ is invertible.
By Definition~\ref{def:TCS}, it is straightforward to verify that

\begin{eqnarray*}
\E{\mS\mS^\top} &=& \begin{bmatrix} \frac{(1-b)\tau_d}{n}\mI_d & 0 \\ 0 & \frac{b\tau_n}{n}\mI_n \end{bmatrix}
\end{eqnarray*}
is invertible and ${\bf Ker}(\E{\mS\mS^\top}) = \{0\}$.
Now we show the Jacobian $DF^\top(x)$ invertible. Let $x = [\alpha; w] \in \R^{n+d}$ with $\alpha \in \R^n$ and $w \in \R^d$. Then $DF(x)$ is written as

\begin{eqnarray} \label{eq:DFalphax}
DF(x)^\top &=&
\begin{bmatrix}
\frac{1}{\lambda n}\mA & - \mI_d \\
\mI_n & \nabla \Phi(w)^\top
\end{bmatrix},
\end{eqnarray}
where $\nabla \Phi(w)^\top = \Diag{\phi_1''(a_1^\top w), \ldots, \phi_n''(a_n^\top w)} \,\mA^\top \in \R^{n \times d}.$
Denote the diagonal matrix $D(w) \eqdef \Diag{\phi_1''(a_1^\top w), \ldots, \phi_n''(a_n^\top w)}$. Since $\phi_i$ is continuously twice differentiable and convex, $\phi_i''(a_i^\top w) \geq 0$ for all $i$. Thus, $D(w) \geq 0$.

Let $(u; v) \in \R^{n+d}$ with $u \in \R^n$ and $v \in \R^d$ such that $DF(x)^\top [u; v] = 0$. We have

\begin{align}
DF(x)^\top \begin{bmatrix}u \\ v\end{bmatrix} = 0 \overset{\eqref{eq:DFalphax}}{\Longleftrightarrow}
\begin{bmatrix}
\frac{1}{\lambda n}\mA & - \mI_d \\
\mI_n & \nabla \Phi(w)^\top
\end{bmatrix} \begin{bmatrix}u \\ v\end{bmatrix} = 0
\Longrightarrow \left(\mI_n + \frac{1}{\lambda n}D(w)\mA^\top \mA\right) u = 0. \label{eq:IDAA}
\end{align}

If $D(w)$ is invertible, \eqref{eq:IDAA} becomes
\begin{eqnarray}
D(w)\left(D(w)^{-1} + \frac{1}{\lambda n} \mA^\top \mA\right) u = 0
&\Longleftrightarrow& \left(D(w)^{-1} + \frac{1}{\lambda n} \mA^\top \mA\right) u = 0. \label{eq:DAA}
\end{eqnarray}
Since $D(w)$ is invertible, i.e.\ $D(w) > 0$, we obtain $D(w)^{-1} > 0$. As $\frac{1}{\lambda n} \mA^\top \mA \geq 0$, we get $D(w)^{-1} + \frac{1}{\lambda n} \mA^\top \mA > 0$, thus invertible. From \eqref{eq:DAA}, we get $u = 0$.

Otherwise, $D(w)$ is not invertible. Without losing generality, we assume that $\phi_1''(a_1^\top w) \geq \phi_2''(a_2^\top w) \geq \cdots \geq \phi_n''(a_n^\top w) = 0$. Let $j$ be the largest index for which $\phi_j''(a_j^\top w) > 0$. If $j$ does not exist, then $D(w) = 0$. From \eqref{eq:IDAA}, we get $u=0$ directly. If $j$ exists, we have $1 \leq j < n$ and
\begin{align}
D(w)\mA^\top \mA \; &= \; \Diag{\phi_1''(a_1^\top w), \ldots, \phi_j''(a_j^\top w), 0, \ldots, 0} \mA^\top \mA \nonumber \\
\; &= \; \begin{bmatrix}
\Diag{\phi_1''(a_1^\top w), \ldots, \phi_j''(a_j^\top w)}\mA_{1:j}^\top \mA_{1:j} & 0 \\
0 & 0 \end{bmatrix}, \label{eq:DAAj}
\end{align}
where $\mA_{1:j} \eqdef [a_1 \ \cdots \ a_j] \in \R^{d \times j}$. Note $u = [u_1; \cdots; u_n] \in \R^n$. Plugging~\eqref{eq:DAAj} into~\eqref{eq:IDAA}, we get
\begin{eqnarray}
& & \left(\mI_n + \frac{1}{\lambda n}
\begin{bmatrix}
\Diag{\phi_1''(a_1^\top w), \ldots, \phi_j''(a_j^\top w)}\mA_{1:j}^\top \mA_{1:j} & 0 \\
0 & 0
\end{bmatrix} \right) u = 0 \nonumber \\
&\Longleftrightarrow& \begin{cases} \left(\mI_j + \frac{1}{\lambda n}\Diag{\phi_1''(a_1^\top w), \ldots, \phi_j''(a_j^\top w)}\mA_{1:j}^\top \mA_{1:j}\right)u_{1:j} = 0 \\ u_{(j+1):n} = 0 \end{cases} \label{eq:IDAAj},
\end{eqnarray}
where $u_{1:j} \eqdef \left[u_1; \cdots; u_j\right] \in \R^j$ and $u_{(j+1):n} \eqdef \left[u_{j+1}; \cdots; u_n\right] \in \R^{n-j}$. From~\eqref{eq:IDAAj}, $u_{(j+1):n} = 0$. Now $\Diag{\phi_1''(a_1^\top w), \ldots, \phi_j''(a_j^\top w)}$ is invertible in the subspace $\R^j$ as every coordinate in the diagonal $\phi_i''(a_i^\top w)$ is strictly positive for all $1 \leq i \leq j$. Similarly, we obtain $u_{1:j} = 0$ from the first equation of \eqref{eq:IDAAj}. Overall we get $u=0$.

Thus, in all cases, $u=0$, then $v = \frac{1}{\lambda n}Au = 0$. We can thus induce that $DF(\alpha; w)^\top$ is invertible for all $\alpha$ and $w$.
\ruiline{Similar to Lemma~\ref{lem:adaptedsketchSNM}, we have~\eqref{eq:adaptedsketch} hold, and by Lemma~\ref{lem:reformulationnew}, we have that~\eqref{eq:DFTinEmHass} holds.}
\end{proof}

\ruiline{From Lemma~\ref{lem:adaptTCS}, we know that for \emph{any} size of the block sketch $\tau_d \in \{1, \cdots, d\}$ and $\tau_n \in \{1, \cdots, n\}$,~\eqref{eq:DFTinEmHass} holds. The corresponding sketch size of $\mS$ is $\tau_d + \tau_n$.}

\subsection{Experiments for TCS method applied for GLM} \label{sec:experiments}

We consider the logistic regression problem with $8$ datasets\footnote{All datasets except for the artificial dataset can be found downloaded on \url{https://www.csie.ntu.edu.tw/~cjlin/libsvmtools/datasets/}.
Some of the datasets can be found originally in~\cite{a9a,covtype,phishing,webspam,UCI}.
} taken from LibSVM~\cite{libsvm}, except for one artificial dataset. Table~\ref{tab:datasets} provides the details of these datasets, including the \emph{condition number} (C.N.) of the model and the smoothness constant $L$ of the model.
C.N.\ of the logistic regression problem is given by
\[\mbox{C.N.} \eqdef \frac{\lambda_{\max}(\mA\mA^\top)}{4n\lambda} + 1,\]
where  $\lambda_{\max}(\cdot)$ is the largest eigenvalue operator.
The smoothness constant $L$ is given by
\[L \eqdef \frac{\lambda_{\max}(\mA\mA^\top)}{4n} + \lambda.\]
As for the logistic regression problem, we consider the loss function $\phi_i$ in \eqref{eq:P} in the form
\[\phi_i(t) = \ln(1+\mathrm{e}^{-y_it})\]
where $y_i$ are the target values for $i = 1, \cdots, n$.

\paragraph{The artificial dataset}
The artificial dataset $\mA^\top \in \R^{n \times d}$ in Table~\ref{tab:datasets} is of size $10000 \times 50$ and generated by a Gaussian distribution whose mean is zero and covariance is  a Toeplitz matrix. Toeplitz matrices are completely determined by their diagonal. We  set the diagonal of our Toeplitz matrix  as
\[[c^0; c^1; \cdots; c^{d-1}] \in \R^d\]
where $c \in \R^+$ is a parameter. We choose $c = 0.9$ (closed to $1$) which results in $\mA$ having highly correlated columns, which in turn makes $\mA$ an  ill-conditioned data set. We set the ground truth coefficients of the model
\[{\bf w} = [(-1)^0 \cdot \mathrm{e}^{-\frac{0}{10}}; \cdots; (-1)^{d-1} \cdot \mathrm{e}^{-\frac{d-1}{10}}] \in \R^d\]
and the target values of the dataset
\[{\bf y} = \text{sgn}\left(\mA^\top {\bf w} + {\bf r}\right) \in \R^n\]
where ${\bf r} \in \R^n$ is the noise generated from a standard normal distribution.

\begin{table}
  \caption{Details of the data sets for binary classification}
  \label{tab:datasets}
  \centering
  \begin{tabular}{ ccccc }
  	 \toprule
     dataset     & dimension ($d$)  & samples ($n$)  & C.N. of the model        & $L$                  \\
     \midrule
     covetype    & $54$             & $581012$       & $7.45 \times 10^{12}$  & $1.28 \times 10^7$          \\
     a9a         & $123$            & $32561$        & $5.12 \times 10^4$     & $1.57$                      \\
     fourclass   & $2$              & $862$          & $4.86 \times 10^6$     & $5.66 \times 10^3$          \\
     artificial  & $50$             & $10000$        & $3.91 \times 10^4$     & $3.91$                      \\
     ijcnn1      & $22$             & $49990$        & $2.88 \times 10^3$     & $5.77 \times 10^{-2}$       \\
     webspam     & $254$            & $350000$       & $7.47 \times 10^4$     & $2.13 \times 10^{-1}$       \\
     epsilon     & $2000$           & $400000$       & $3.51 \times 10^4$     & $8.76 \times 10^{-2}$       \\
     phishing    & $68$             & $11055$        & $1.04 \times 10^3$     & $9.40 \times 10^{-2}$       \\
     \bottomrule
  \end{tabular}
\end{table}

We compare the TCS method with SAG~\cite{SAG}, SVRG~\cite{Johnson2013}, dfSDCA~\cite{Shalev-Shwartz2015} and Quartz~\cite{Qu2015b}. All experiments were initialized at $w^0 = 0 \in \R^d$ (and/or $\alpha^0 = 0 \in \R^n$ for TCS/dfSDCA methods)
 and were performed in Python 3.7.3 on a laptop with an Intel Core i9-9980HK CPU and 32 Gigabyte of DDR4 RAM running OSX 10.14.5.
For all methods, we used the stepsize that was shown to work well in practice. For instance, the common rule of thumb for SAG and SVRG is to use a stepsize $\frac{1}{L}$, where $L$ is the smoothness constant. This rule of thumb stepsize is not supported by theory. Indeed for SAG, the theoretical stepsize is $\frac{1}{16L}$ and it should be even smaller for SVRG depending on the condition number. For dfSDCA and Quartz’s, we used the stepsize suggested in the experiments in~\cite{Shalev-Shwartz2015} and~\cite{Qu2015b} respectively.

For TCS, we used two types of stepsize, related to the C.N.\ of the model. If the condition number is big (Fig.~\ref{figure} top row), we used $\gamma = 1$ except for a9a with $\gamma = 1.5$. If the condition number is small (Fig.~\ref{figure} bottom row), we used $\gamma = 1.8$.
We also set the Bernoulli parameter $b$ (probability of the coin toss) depending on the size of the dataset (see Table~\ref{tab:parameters} in Sec.~\ref{sec:detailexp}), and $\tau_d = d$. We tested three different sketch sizes $\tau_n = 50, 150, 300$. More details of the parameter settings are presented in Sec.~\ref{sec:detailexp}.

We used $\lambda = \frac{1}{n}$ regularization parameter, evaluated each method $10$ times and stopped once the gradient norm\footnote{\xline{We evaluated the true gradient norm every $1000$ iterations. We also paused the timing when computing the performance evaluation of the gradient norm.}} was below $10^{-5}$ or some maximum time had been reached. In Fig.~\ref{figure}, we plotted the central tendency as a solid line and all other executions as a shaded region for the wall-clock time vs gradient norm.

From Fig.~\ref{figure}, TCS outperforms all other methods on ill-conditioned problems (Fig.~\ref{figure} top row), but not always the case on well-conditioned problems (Fig.~\ref{figure} bottom row). This is because in ill-conditioned problems, the curvature of the optimization function is not uniform over directions and varies in the input space. Second-order methods effectively exploit information of the local curvature to converge faster than first-order methods.
To further illustrate the performance of TCS on ill-conditioned problems, we compared the performance of TCS on the artificial dataset

in the top right of Fig.~\ref{figure}.
Note as well that for reaching an approximate solution at early stage (i.e.\ $tol = 10^{-3}, 10^{-4}$), TCS is very competitive on all problems. TCS also has the smallest variance compared to the first-order methods based on eye-balling the shaded error bars in Fig.~\ref{figure}, especially compared to SVRG. Among the three tested sketch sizes, $150$ performed the best except on the \textit{epsilon} dataset.

\begin{figure}
\centering
\begin{tabular}{cccc}
\includegraphics[width=.21\linewidth]{img/covtype_gnorm.pdf}&
\includegraphics[width=.21\linewidth]{img/a9a_gnorm.pdf}&
\includegraphics[width=.21\linewidth]{img/fourclass_gnorm.pdf}&
\includegraphics[width=.21\linewidth]{img/art_gnorm.pdf} \\
\includegraphics[width=.21\linewidth]{img/ijcnn1_gnorm.pdf}&
\includegraphics[width=.21\linewidth]{img/webspam_gnorm.pdf}&
\includegraphics[width=.21\linewidth]{img/epsilon_normalized_gnorm.pdf}&
\includegraphics[width=.21\linewidth]{img/phishing_gnorm.pdf}
\end{tabular}
\caption{Experiments for TCS method applied for generalized linear model.}
\label{figure}
\end{figure}

\section{Conclusion and future work}

We introduced the \SNR method, for which we provided strong convergence guarantees. We also developed several promising applications of \SNR to show that \SNR is very flexible and tested one of these specialized variants for training GLMs. \SNR is flexible by the fact that its primitive goal is to solve efficiently nonlinear equations. Since there are many ways to re-write an optimization problem as nonlinear equations, each re-write leads to a distinct method, thus leads to a specific implementation in practice (e.g.\ \SNM, TCS methods) when using \SNR. Besides, the convergence theories presented in Sec.~\ref{sec:convergence} guarantee a large variety of choices for the sketch. This flexibility allows us to discover many applications of \SNR and their induced consequences, especially providing new global convergence theories. As such, we believe that \SNR and its global convergence theory will open the way to designing and analyzing a host of new stochastic second order methods. Further venues of investigation include exploring the use of adaptive norms for projections and leveraging efficient sketches (e.g.\ the fast Johnson-Lindenstrauss sketch~\cite{Pilanci2015a}, sketches with determinantal sampling~\cite{mutny2020convergence}) to design even faster variants of \SNR or cover other stochastic second-order methods. Since \SNR can be seen as \SGD, it might be possible to design and develop efficient accelerated \SNR or \SNR with momentum methods. On the experimental side, it would be interesting to apply our method to the training of deep neural networks.