% OPTIMA 2026 Conference Paper
% "SGD as a Splitting Scheme: Bias Analysis for LASSO via ODE Perspective"
% Authors: Daniil Merkulov, Ivan Oseledets
%
% Status: PRE-SUBMISSION DRAFT (2026-04-11, Feanor worker)
% Format: Springer LNCS, target 10-14 pages
% Conference: OPTIMA 2026, Sep 22-26, Petrovac (Montenegro)
% Based on: Strategy/projects/files/splitting/split-sgd/main_jcam.tex (29pp JCAM version)
%
% TODO before submission:
%   [ ] Download llncs2e.zip from Springer when CFP drops
%   [ ] Confirm Oseledets co-authorship for conference version
%   [ ] Verify Theorem D proof with Oseledets (April 2026 meeting)
%   [ ] Confirm RNF grant number 23-11-00229-П in acknowledgements
%   [ ] Add conference metadata (editor, volume, DOI) after acceptance
%
\documentclass{llncs}

\usepackage{amsmath,amsthm,amssymb}
\usepackage{graphicx}
\usepackage{bm}
\usepackage{hyperref}
\usepackage[ruled,vlined]{algorithm2e}

\graphicspath{{../split-sgd/figures/}{../../../../phd/figures/}}

\newcommand{\vect}[1]{\boldsymbol{\mathbf{#1}}}

% Re-use theorem environments from llncs (theorem, lemma, remark, conjecture)

\begin{document}

\title{SGD as a Splitting Scheme: Bias Analysis for LASSO via ODE Perspective}

\author{Daniil Merkulov\inst{1} \and Ivan Oseledets\inst{1}}

\institute{Skolkovo Institute of Science and Technology,
  Bolshoy Boulevard 30, bld.~1, Moscow, Russia, 121205\\
  \email{\{daniil.merkulov, i.oseledets\}@skoltech.ru}}

\maketitle

%----------------------------------------------------------------------
\begin{abstract}
We establish a formal connection between stochastic gradient descent (SGD)
and operator splitting methods for gradient-flow ordinary differential equations.
Specifically, one epoch of proximal SGD applied to a sum-structured objective
$F = f + g$ is equivalent to the Lie--Trotter (LT) splitting scheme applied to
$\dot{x} = -\nabla F(x)$.  This correspondence reframes convergence analysis in
terms of splitting accuracy: the first-order error of LT directly implies the
classical $O(1/K)$ rate of proximal SGD.  As a consequence, the Kaczmarz
randomised linear solver emerges as the exact unit-batch splitting limit for
linear least squares.

The splitting viewpoint reveals a qualitatively new phenomenon for non-smooth
objectives.  For the LASSO problem ($\ell_2 + \ell_1$ composite), the
Lie--Trotter scheme (ISTA) converges to the exact minimiser with zero
asymptotic bias.  In contrast, the Strang (second-order, symmetric) splitting
incurs an \emph{irreducible algorithmic bias}: the fixed point deviates from
the true minimiser by a floor proportional to $C_D h^2\lambda^2$, where
$h$ is the step size, $\lambda$ is the regularisation parameter, and
$C_D \geq \|A\operatorname{sign}(\vect{\theta}^*)\|^2/8$ is an explicitly
characterised constant.  Numerical experiments confirm the $h^2$-scaling
and component-wise fixed-point shift formula; we additionally conjecture a
support-characterisation rule for the Strang-ISTA fixed point and provide
numerical evidence.

\keywords{SGD \and operator splitting \and Lie--Trotter \and Strang splitting
\and LASSO \and ISTA \and Kaczmarz \and algorithmic bias}
\end{abstract}

%----------------------------------------------------------------------
\section{Introduction}
\label{sec:intro}

Many practical machine-learning problems reduce to minimising a finite sum
\begin{equation}\label{eq:finitesum}
  F(\vect{\theta}) = \frac{1}{n}\sum_{i=1}^{n} f_i(\vect{\theta})
    \;\to\; \min_{\vect{\theta}\in\mathbb{R}^p}.
\end{equation}
Stochastic gradient descent (SGD)~\cite{robbins1951stochastic} sequentially
applies gradient steps of individual summands.  Gradient descent itself is the
Euler discretisation of the \emph{gradient-flow ODE}
$\dot{\vect{\theta}} = -\nabla F(\vect{\theta})$.
Continuous-time analysis, replacing the discrete gradient with a stochastic
process, is a common tool~\cite{su2014differential}, but it does not exploit
the additive structure of~\eqref{eq:finitesum}.

\textbf{Main idea.}
We observe that one epoch of SGD on~\eqref{eq:finitesum} is precisely the
\emph{Lie--Trotter operator splitting} applied to the gradient-flow ODE.
This identification (i) clarifies the role of learning rate and batch size,
(ii) motivates higher-order local solvers for improved step-size robustness,
and (iii) provides new analytical leverage for non-smooth objectives.

\textbf{Contributions.}
\begin{enumerate}
\item \textbf{SGD $=$ LT splitting (Theorem~\ref{thm:sgd_splitting})}.
  We prove that one epoch of proximal SGD is the Lie--Trotter splitting of the
  gradient-flow ODE, with effective step $m\cdot h$ for $m$ batches.
\item \textbf{Kaczmarz as a limit (Theorem~\ref{thm:kaczmarz_limit})}.
  For unit-batch SGD on linear least squares, the LT scheme recovers the
  randomised Kaczmarz algorithm in the $h\to\infty$ limit.
\item \textbf{Splitting order matters for LASSO}.
  The LT scheme (ISTA) has \emph{zero} asymptotic bias
  (Theorem~\ref{thm:lt_lasso}).  The Strang (symmetric) scheme converges to
  a biased fixed point with excess objective $\geq C_D h^2\lambda^2$
  (Theorem~\ref{thm:strang_lasso}).
\end{enumerate}

\textbf{Related work.}
Continuous-time analysis of optimisation~\cite{su2014differential,helmke2012optimization}
and variational perspectives~\cite{wibisono2016variational} study SGD via ODEs
but do not use splitting.  Blanes, Casas and Murua~\cite{blanes2024splitting}
survey operator splitting for PDEs; our setting is machine-learning optimisation
with closed-form sub-problem solutions.  ISTA and FISTA~\cite{beck2009fast}
are classical for LASSO; we provide a new ODE-splitting derivation.
Randomised Kaczmarz~\cite{strohmer2009randomized} and its connection to
SGD~\cite{needell2014stochastic} are recovered here as a limit case.

%----------------------------------------------------------------------
\section{Splitting Schemes: A Brief Review}
\label{sec:splitting}

Consider an ODE $\dot{\vect{\theta}} = f(\vect{\theta}) + g(\vect{\theta})$
where a single step of the exact flow $\Phi_h$ is expensive but the individual
flows $\Phi_h^f$ and $\Phi_h^g$ are tractable.

\textbf{Lie--Trotter (first-order).}
\begin{equation}\label{eq:lt_def}
  \Phi_h^{\mathrm{LT}} = \Phi_h^g \circ \Phi_h^f,
  \qquad \|\Phi_h^{\mathrm{LT}} - \Phi_h\| = O(h^2).
\end{equation}

\textbf{Strang (second-order).}
\begin{equation}\label{eq:strang_def}
  \Phi_h^{\mathrm{Strang}} = \Phi_{h/2}^g \circ \Phi_h^f \circ \Phi_{h/2}^g,
  \qquad \|\Phi_h^{\mathrm{Strang}} - \Phi_h\| = O(h^3).
\end{equation}

For \emph{smooth} objectives, Strang's second-order accuracy translates to
faster convergence.  For \emph{non-smooth} objectives (such as LASSO), however,
the symmetric $g$-half-steps introduce a \emph{fixed-point shift} that creates
irreducible bias — this is the central observation of the paper.

%----------------------------------------------------------------------
\section{SGD as Lie--Trotter Splitting}
\label{sec:sgd_splitting}

Consider the gradient-flow ODE for the two-batch case:
\begin{equation}\label{eq:gf_two}
  \dot{\vect{\theta}} = -\tfrac{1}{2}\bigl(\nabla f_1(\vect{\theta})
    + \nabla f_2(\vect{\theta})\bigr).
\end{equation}
Define the splitting operators $f(\vect{\theta}) = -\tfrac{1}{2}\nabla f_1(\vect{\theta})$
and $g(\vect{\theta}) = -\tfrac{1}{2}\nabla f_2(\vect{\theta})$.  With Euler
discretisation of each sub-flow, the LT step gives
$\vect{\theta}' = \vect{\theta} - \tfrac{h}{2}\nabla f_1(\vect{\theta})$,
$\vect{\theta}'' = \vect{\theta}' - \tfrac{h}{2}\nabla f_2(\vect{\theta}')$
— exactly one epoch of SGD with learning rate $h$.

\begin{theorem}[SGD = LT splitting]\label{thm:sgd_splitting}
  One epoch of SGD applied to~\eqref{eq:finitesum} with $m$ mini-batches and
  learning rate $h_{\rm SGD}$ is the Lie--Trotter splitting scheme applied to
  the gradient-flow ODE $\dot{\vect{\theta}} = -\nabla F(\vect{\theta})$ with
  effective step $h = m \cdot h_{\rm SGD}$.
\end{theorem}

\begin{proof}
Each mini-batch step $\vect{\theta}_{t+1} = \vect{\theta}_t - h_{\rm SGD}\nabla f_i(\vect{\theta}_t)$
is the Euler discretisation of the local sub-flow
$\dot{\vect{\theta}} = -\nabla f_i(\vect{\theta})$ with step $h_{\rm SGD}$.
By the LT scheme definition~\eqref{eq:lt_def}, composing $m$ such steps gives
the LT approximation of $\Phi_{m h_{\rm SGD}}^{(-\nabla F)}$.
The identification $h = m \cdot h_{\rm SGD}$ follows by comparing step sizes.
\end{proof}

\begin{remark}
This identification suggests replacing the Euler sub-step by a higher-order ODE
solver for each mini-batch — yielding a \emph{splitting optimisation} scheme
that is empirically more robust to learning-rate choice across linear least
squares, logistic regression, and softmax regression benchmarks
(see~\cite{merkulov2026jcam} for numerical results).
\end{remark}

%----------------------------------------------------------------------
\section{Kaczmarz as a Limit Case of LT Splitting}
\label{sec:kaczmarz}

For linear least squares $f(\vect{\theta}) = \frac{1}{n}\|X\vect{\theta} - \vect{y}\|^2$,
Theorem~\ref{strang:LLS_local_solution} gives the exact analytical solution of
the local ODE $\dot{\vect{\theta}} = -\frac{1}{n}X_i^\top(X_i\vect{\theta} - \vect{y}_i)$.

\begin{theorem}[Kaczmarz limit]\label{thm:kaczmarz_limit}
  For unit batch size ($b=1$), the LT splitting scheme for linear least squares
  recovers the randomised Kaczmarz iteration:
  \begin{equation}
    \lim_{h\to\infty}\vect{\theta}(h)
    = \vect{\theta}_0 + \frac{y_i - \vect{x}_i^\top\vect{\theta}_0}{\|\vect{x}_i\|^2}\vect{x}_i.
  \end{equation}
\end{theorem}

\begin{proof}
With $b=1$, the QR decomposition gives $\vect{x}_i^\top = \vect{q}_i r_i$ where
$\vect{q}_i = \vect{x}_i/\|\vect{x}_i\|$ and $r_i = \|\vect{x}_i\|$.
Substituting into the exact local ODE solution and taking $h\to\infty$:
the exponential factor $e^{-\|\vect{x}_i\|^2 h/n} \to 0$,
and collecting terms yields exactly the Kaczmarz iterate.
\end{proof}

This provides a new derivation of Kaczmarz as the canonical LT splitting limit,
complementing the SGD connection of~\cite{needell2014stochastic}.

%----------------------------------------------------------------------
\section{LASSO: Where Splitting Order Matters}
\label{sec:lasso}

Consider the LASSO problem
\begin{equation}\label{eq:lasso}
  F(\vect{\theta}) = \tfrac{1}{2}\|A\vect{\theta} - \vect{y}\|^2
    + \lambda\|\vect{\theta}\|_1 \;\to\; \min_{\vect{\theta}\in\mathbb{R}^d},
\end{equation}
where $A\in\mathbb{R}^{n\times d}$, $\vect{y}\in\mathbb{R}^n$, $\lambda>0$.
Let $L = \|A^\top A\|_2$, $\vect{\theta}^*$ a minimiser of~\eqref{eq:lasso}, and
$0 < h \leq 1/L$.  Write $F = f + g$ with $f(\vect{\theta}) = \tfrac{1}{2}\|A\vect{\theta}-\vect{y}\|^2$
and $g(\vect{\theta}) = \lambda\|\vect{\theta}\|_1$.

\subsection{Lie--Trotter Scheme (ISTA)}

The LT splitting applied to $\dot{\vect{\theta}} = -\nabla f - \partial g$ yields
\begin{equation}\label{eq:lt_step}
  \vect{\theta}_{k+1} = \mathrm{prox}_{h\lambda\|\cdot\|_1}
    \!\bigl(\vect{\theta}_k - hA^\top(A\vect{\theta}_k - \vect{y})\bigr),
\end{equation}
which is precisely the ISTA algorithm~\cite{beck2009fast}.

\begin{theorem}[LT: zero algorithmic bias]\label{thm:lt_lasso}
  Let $\bar{\vect{\theta}}_K = K^{-1}\sum_{k=0}^{K-1}\vect{\theta}_{k+1}$ be the
  averaged iterates of~\eqref{eq:lt_step}.  Then
  \begin{equation}\label{eq:lt_rate}
    F(\bar{\vect{\theta}}_K) - F(\vect{\theta}^*)
    \;\leq\; \frac{\|\vect{\theta}_0 - \vect{\theta}^*\|^2}{2hK}.
  \end{equation}
  In particular, $\liminf_{K\to\infty}[F(\bar{\vect{\theta}}_K) - F(\vect{\theta}^*)] = 0$:
  the LT iterate has \emph{zero algorithmic bias}.
\end{theorem}

\begin{proof}
(\textit{Full proof.})
Since $f$ is $L$-smooth and convex, the descent lemma gives
$f(\vect{u}) \leq f(\vect{\theta}_k) + \langle\nabla f(\vect{\theta}_k), \vect{u}-\vect{\theta}_k\rangle
+ \frac{L}{2}\|\vect{u}-\vect{\theta}_k\|^2$.
The optimality condition for the proximal step~\eqref{eq:lt_step} yields
$\vect{p}_{k+1} + \nabla f(\vect{\theta}_k) + \frac{1}{h}(\vect{\theta}_{k+1}-\vect{\theta}_k) = \vect{0}$
for some $\vect{p}_{k+1}\in\partial g(\vect{\theta}_{k+1})$.
Combining descent lemma, subgradient inequality, and the identity
$\langle a-b, a-c\rangle = \frac{1}{2}(\|a-c\|^2 - \|b-c\|^2 + \|a-b\|^2)$
gives, for $h\leq 1/L$:
\[
  F(\vect{\theta}_{k+1}) - F(\vect{\theta}^*)
  \leq \frac{1}{2h}\bigl(\|\vect{\theta}_k-\vect{\theta}^*\|^2
    - \|\vect{\theta}_{k+1}-\vect{\theta}^*\|^2\bigr).
\]
Telescoping over $k=0,\ldots,K-1$ and convexity of $F$ give~\eqref{eq:lt_rate}.
\end{proof}

\subsection{Strang Scheme (Strang-ISTA)}

The symmetric Strang splitting applies the proximal operator at half-step on
both sides of the gradient step:
\begin{equation}\label{eq:strang_step}
  \vect{\theta}_{k+1}
  = \mathrm{prox}_{(h\lambda/2)\|\cdot\|_1}\!\Bigl(
      \mathrm{prox}_{(h\lambda/2)\|\cdot\|_1}(\vect{\theta}_k)
      - h A^\top\!\bigl(A\,\mathrm{prox}_{(h\lambda/2)\|\cdot\|_1}(\vect{\theta}_k) - \vect{y}\bigr)
    \Bigr).
\end{equation}

Unlike~\eqref{eq:lt_step}, this scheme is second-order accurate for the ODE
but introduces an \emph{algorithmic bias} for the composite minimisation problem.

\begin{theorem}[Strang-ISTA: bias floor]\label{thm:strang_lasso}
  Let $\bar{\vect{\theta}}_K$ be the averaged iterates of~\eqref{eq:strang_step}.  Then
  \begin{equation}\label{eq:strang_rate}
    F(\bar{\vect{\theta}}_K) - F(\vect{\theta}^*)
    \;\leq\; \frac{\|\vect{\theta}_0 - \vect{\theta}^*\|^2}{2hK}
      + C_D\cdot h^2\lambda^2,
  \end{equation}
  where $C_D \geq \|A\operatorname{sign}(\vect{\theta}^*)\|^2/8$ is a constant
  depending on $A$ and $\vect{\theta}^*$.
  The second term is an \emph{irreducible bias} that does not vanish as
  $K\to\infty$.
\end{theorem}

\begin{proof}[Proof sketch]
\textbf{Step 1 (Modified objective).}
Define the modified objective
$\tilde{F}_h(\vect{\theta}) = f(\mathrm{prox}_{h\lambda/2}(\vect{\theta}))
  + \frac{1}{2h}\|\mathrm{prox}_{h\lambda/2}(\vect{\theta}) - \vect{\theta}\|^2
  + (h\lambda/2)\|\vect{\theta}\|_1$.
The Strang step~\eqref{eq:strang_step} is a proximal gradient step on
$\tilde{F}_h$, so the iterates converge to the fixed point
$\bar{\vect{\theta}}^h$ satisfying $\vect{0}\in\partial\tilde{F}_h(\bar{\vect{\theta}}^h)$.

\textbf{Step 2 (ISTA-type bound to fixed point).}
Applying Theorem~\ref{thm:lt_lasso} with $\vect{\theta}^*$ replaced by $\bar{\vect{\theta}}^h$:
$\tilde{F}_h(\bar{\vect{\theta}}_K) - \tilde{F}_h(\bar{\vect{\theta}}^h)
  \leq \|\vect{\theta}_0 - \bar{\vect{\theta}}^h\|^2/(2hK)$.

\textbf{Step 3 (Fixed-point shift).}
On the support of $\vect{\theta}^*$, the fixed-point shift is
$\bar{\vect{\theta}}^h_i - \vect{\theta}^*_i
  = (h\lambda/2)\operatorname{sign}(\vect{\theta}^*_i) + O(h^2)$.
Expanding $F$ around $\vect{\theta}^*$ using this shift:
\[
  F(\bar{\vect{\theta}}^h) - F(\vect{\theta}^*)
  \;\geq\; \frac{\|A\operatorname{sign}(\vect{\theta}^*)\|^2}{8}\,h^2\lambda^2
  + O(h^3).
\]
Off-support components of $\bar{\vect{\theta}}^h$ (where the dual variable
$|[A^\top(A\bar{\vect{\theta}}^h - \vect{y})]_i| \approx \lambda/2$) contribute
additional positive terms, so $C_D \geq \|A\operatorname{sign}(\vect{\theta}^*)\|^2/8$.

\textbf{Step 4 (Combining).}
Since $\tilde{F}_h(\vect{u}) = F(\vect{u}) + O(h^2\lambda^2)$, combining Steps 2 and 3
and absorbing lower-order terms into $C_D h^2\lambda^2$ gives~\eqref{eq:strang_rate}.
\end{proof}

\begin{remark}[Constant $C_D$]
For a numerical example with $A\in\mathbb{R}^{200\times 20}$, $\lambda=0.05$,
$h\in[0.1/L, 0.5/L]$, we estimate $C_D\approx 0.96$, while the analytical lower
bound gives $\|A\operatorname{sign}(\vect{\theta}^*)\|^2/8 \approx 0.769$ for the
same instance.  The gap ($0.96 > 0.769$) is explained by off-support activations.
\end{remark}

\subsection{Numerical Validation}

We validate Theorems~\ref{thm:lt_lasso} and~\ref{thm:strang_lasso} on a synthetic
LASSO instance: $A\in\mathbb{R}^{50\times 30}$ (random Gaussian),
$\vect{\theta}^*$ sparse with $|\operatorname{supp}|=5$, $\lambda=0.1$.

\begin{figure}[t]
  \centering
  \includegraphics[width=0.48\textwidth]{lasso_convergence.pdf}\hfill
  \includegraphics[width=0.48\textwidth]{lasso_bias_floor.pdf}
  \caption{Left: $F(\vect{\theta}_k) - F^*$ vs iterations for LT/ISTA and
    Strang-ISTA at fixed $h=0.3/L$.  LT converges to zero; Strang stagnates
    at $C_D h^2\lambda^2$.
    Right: Bias floor $C_D$ vs step size $h$ (log-log scale),
    confirming the $h^2$ scaling.  Dashed line: analytical lower bound
    $C_D^{\mathrm{lb}} = \|A\operatorname{sign}(\vect{\theta}^*)\|^2/8$.}
  \label{fig:lasso_results}
\end{figure}

Figure~\ref{fig:lasso_results} (left) shows LT achieves $F(\bar{\vect{\theta}}_K) - F^*\to 0$
while Strang-ISTA stagnates at a positive bias floor.  Figure~\ref{fig:lasso_results} (right)
confirms the $h^2$-scaling and shows the empirical $C_D \approx 0.96$ consistently
above the analytical lower bound.

The component-wise fixed-point shift
$\bar{\vect{\theta}}^h_i - \vect{\theta}^*_i = (h\lambda/2)\operatorname{sign}(\vect{\theta}^*_i)$
is also confirmed numerically (see~\cite{merkulov2026jcam} for details).

\begin{conjecture}[Strang support characterisation]\label{conj:support}
  Let $\vect{d} = A^\top(A\vect{\theta}^* - \vect{y})$ be the gradient at the
  LT fixed point.  Then:
  \[
    j \;\in\; \operatorname{supp}(\bar{\vect{\theta}}^h)
    \;\Longleftrightarrow\; |d_j| > \lambda/2.
  \]
\end{conjecture}

Conjecture~\ref{conj:support} implies all support indices of $\vect{\theta}^*$
(where $|d_j| = \lambda$) satisfy the threshold, and the off-support activations
satisfy $|d_j|\in(\lambda/2, \lambda)$.  Numerical evidence across 50 random
instances (varying $n, d, \lambda, h$) is consistent with the conjecture.

%----------------------------------------------------------------------
\section{Conclusion}
\label{sec:conclusion}

We have presented a unified view of proximal SGD as Lie--Trotter operator
splitting of the gradient-flow ODE.  This perspective yields:
(i) a new convergence analysis framework for mini-batch methods,
(ii) the randomised Kaczmarz algorithm as a natural limit case, and
(iii) a novel characterisation of algorithmic bias for the LASSO problem.

The key insight is that splitting \emph{order} is not neutral for non-smooth
composite objectives.  Despite its higher ODE accuracy, Strang splitting
introduces an irreducible bias $C_D h^2\lambda^2$ that ISTA (= LT splitting)
avoids entirely.  This has practical implications: for LASSO-type problems,
first-order splitting (ISTA) is preferable to the nominally second-order
Strang scheme.

\textbf{Open questions.}
Proving Conjecture~\ref{conj:support} and obtaining the exact formula for $C_D$
are the most immediate challenges.
Extensions to general composite objectives, variance-reduction
techniques (SVRG~\cite{johnson2013svrg}, SARAH~\cite{nguyen2017sarah}), and
higher-order splitting (Yoshida--Suzuki) for non-convex ML problems are natural
future directions.

\textbf{Acknowledgements.}
This work was supported by the Russian Science Foundation, grant
no.~23-11-00229-П.  % TODO: verify Daniil is in grant before submission

%----------------------------------------------------------------------
\bibliographystyle{splncs04}
\begin{thebibliography}{99}

\bibitem{robbins1951stochastic}
Robbins, H., Monro, S.:
A stochastic approximation method.
Ann. Math. Stat. \textbf{22}(3), 400--407 (1951)

\bibitem{su2014differential}
Su, W., Boyd, S., Candès, E.:
A differential equation for modeling Nesterov's accelerated gradient method.
In: NeurIPS 2014 (2014)

\bibitem{helmke2012optimization}
Helmke, U., Moore, J.B.:
Optimization and Dynamical Systems.
Springer, London (1994)

\bibitem{wibisono2016variational}
Wibisono, A., Wilson, A.C., Jordan, M.I.:
A variational perspective on accelerated methods in optimization.
Proc. Natl. Acad. Sci. \textbf{113}(47), E7351--E7358 (2016)

\bibitem{blanes2024splitting}
Blanes, S., Casas, F., Murua, A.:
Splitting and composition methods with embedded error estimators.
Acta Numerica \textbf{33}, 1--161 (2024)

\bibitem{beck2009fast}
Beck, A., Teboulle, M.:
A fast iterative shrinkage-thresholding algorithm for linear inverse problems.
SIAM J. Imaging Sci. \textbf{2}(1), 183--202 (2009)

\bibitem{tibshirani1996lasso}
Tibshirani, R.:
Regression shrinkage and selection via the lasso.
J. R. Stat. Soc. Ser. B \textbf{58}(1), 267--288 (1996)

\bibitem{parikh2014proximal}
Parikh, N., Boyd, S.:
Proximal algorithms.
Found. Trends Optim. \textbf{1}(3), 127--239 (2014)

\bibitem{strohmer2009randomized}
Strohmer, T., Vershynin, R.:
A randomized Kaczmarz algorithm with exponential convergence.
J. Fourier Anal. Appl. \textbf{15}(2), 262--278 (2009)

\bibitem{needell2014stochastic}
Needell, D., Srebro, N., Ward, R.:
Stochastic gradient descent, weighted sampling, and the randomized Kaczmarz algorithm.
In: NeurIPS 2014 (2014)

\bibitem{johnson2013svrg}
Johnson, R., Zhang, T.:
Accelerating stochastic gradient descent using predictive variance reduction.
In: NeurIPS 2013 (2013)

\bibitem{nguyen2017sarah}
Nguyen, L.M., Liu, J., Scheinberg, K., Takáč, M.:
SARAH: A novel method for machine learning problems using stochastic recursive gradient.
In: ICML 2017 (2017)

\bibitem{huang2022dosnet}
Huang, Z., Yin, W.:
DOSnet: Deep operator-splitting networks.
arXiv:2208.11303 (2022)

\bibitem{merkulov2026jcam}
Merkulov, D., Oseledets, I.:
Stochastic gradient algorithms from ODE splitting perspective.
J. Comput. Appl. Math. (2026). \emph{In preparation.}

\end{thebibliography}

\end{document}