\documentclass[preprint,12pt]{elsarticle}

\usepackage{microtype}
\usepackage{graphicx}
\usepackage{subcaption}
\usepackage{booktabs}
\usepackage{hyperref}
\graphicspath{{../figures/}{../../../../phd/figures/}}
\usepackage{amsmath,amsthm,amssymb}
\usepackage{thmtools, thm-restate}
\newtheorem{theorem}{Theorem}
\newtheorem{lemma}{Lemma}
\newtheorem{remark}{Remark}
\newtheorem{corollary}{Corollary}
\newtheorem{conjecture}{Conjecture}
\usepackage{bm}
\newcommand{\vect}[1]{\boldsymbol{\mathbf{#1}}}

\newcommand{\theHalgorithm}{\arabic{algorithm}}
\usepackage[ruled,vlined]{algorithm2e}

\journal{Journal of Computational and Applied Mathematics}

\begin{document}

\begin{frontmatter}

\title{Stochastic gradient algorithms from ODE splitting perspective}

\author[label1]{Daniil Merkulov}
\ead{daniil.merkulov@skolkovotech.ru}

\author[label1]{Ivan Oseledets}
\ead{i.oseledets@skoltech.ru}

\address[label1]{Center for Computational and Data-Intensive Science and Engineering,
Skolkovo Institute of Science and Technology,
Bolshoy Boulevard 30, bld. 1, Moscow, Russia, 121205}

\begin{abstract}
We present a unified view of stochastic gradient methods through the lens of operator splitting for gradient-flow ODEs.
We show that one epoch of SGD is precisely a first-order Lie--Trotter splitting of the gradient-flow ODE, derive a new upper bound on the global splitting error, and establish that the Kaczmarz method is the exact limiting case of unit-batch SGD for the linear least-squares problem.
Replacing the Euler step by a higher-order ODE solver at each mini-batch yields a \emph{splitting optimisation} scheme that is markedly more robust to the learning-rate choice across linear least squares, logistic regression, and softmax regression benchmarks.
Beyond smooth objectives, we analyse the LASSO ($\ell_2 + \ell_1$) problem under both the first-order Lie--Trotter (ISTA) and second-order Strang splitting (Strang-ISTA):
the first-order scheme converges to the exact sparse minimiser with zero asymptotic bias,
whereas the Strang scheme converges to a biased fixed point whose excess objective scales as $C_D h^2\lambda^2$,
with $C_D \geq \|A\operatorname{sign}(\vect{\theta}^*)\|^2/8$.
Numerical experiments confirm the $h^2$-scaling of the bias floor and the component-wise fixed-point shift formula; we additionally conjecture a support-characterisation criterion for Strang-ISTA fixed points ($j \in \operatorname{supp}(\bar{\vect{\theta}}^h) \Leftrightarrow |d_j| > \lambda/2$) and provide numerical evidence in its favour.
\end{abstract}

\begin{keyword}
stochastic gradient descent \sep operator splitting \sep ODE methods \sep
LASSO \sep ISTA \sep Kaczmarz method \sep numerical optimization
\end{keyword}

\begin{highlights}
\item SGD reinterpreted as Lie--Trotter operator splitting of gradient-flow ODE
\item Higher-order local ODE solvers yield step-size-robust splitting optimisation
\item Kaczmarz method is the exact unit-batch splitting limit for linear least squares
\item Lie--Trotter ISTA: zero asymptotic bias; Strang-ISTA has $O(h^2\lambda^2)$ bias floor
\item Strang-ISTA support conjectured: $j \in \operatorname{supp}(\bar{\vect{\theta}}^h) \Leftrightarrow |d_j|>\lambda/2$; numerically confirmed
\end{highlights}

\end{frontmatter}

\section{Introduction}
A lot of practical problems arising in machine learning require minimization of a finite sample average which can be written in the form
\begin{equation}\label{strang:finitesum}
    f(\vect{\vect{\theta}}) = \frac{1}{n} \sum_{i=1}^n f_i(\vect{\vect{\theta}}) \rightarrow \min_{\vect{\vect{\theta}} \in \mathbb{R}^p},
\end{equation}
where the sum goes over the \emph{minibatches} of the original dataset. Vanilla stochastic gradient descent
(SGD) method \cite{robbins1951stochastic} consists sequential steps in the direction of the gradient of $f_i(\vect{\theta})$, where $i$ is to be chosen randomly from $1$ to $n$ without replacement.
\begin{equation}
\vect{\theta}_{k+1} = \vect{\theta}_{k} - h_{k} \nabla f_i.
\end{equation}
Gradient descent method \cite{cauchy1847methode} can be considered as an Euler discretization of the ordinary
differential equation (ODE) of the form of the gradient flow
\begin{equation}\label{strang:euler}
    \frac{d \vect{\theta}}{d t} = -\nabla f(\vect{\theta}).
\end{equation}
In continuous time, SGD if often analyzed by introducing noise into the right-hand side of \eqref{strang:euler}. However, for a real dataset, the distribution of the noise obtained by replacing the full gradient by its minibatch variant is not known and can be different for different problems. Instead, we propose a new view on the SGD as a \emph{first-order splitting scheme} for \eqref{strang:euler}, thus shedding a new light on SGD-type algorithms. This representation allows using more efficient local problem solvers for the approximation of the full gradient flow.

\textbf{Contributions}
\begin{itemize}
    \item We show that vanilla SGD is a first-order Lie--Trotter splitting scheme for the full gradient-flow ODE, clarifying the role of learning rate and batch size in the continuous-time limit.
    \item We propose a splitting optimisation scheme that replaces the Euler step by a higher-order ODE solver at each mini-batch and show empirically that it is a stepsize-robust alternative to SGD for practical ML problems.
    \item We establish that the Kaczmarz (ART) method is the exact limiting case of the unit-batch splitting scheme for the linear least-squares problem.
    \item We analyse the LASSO ($\ell_2 + \ell_1$) problem under both Lie--Trotter (ISTA) and Strang splitting (Strang-ISTA): the first-order scheme achieves zero asymptotic bias, while the Strang scheme incurs a bias floor scaling as $C_D h^2\lambda^2$. We characterise the component-wise fixed-point shift analytically, conjecture a support-characterisation rule for Strang-ISTA fixed points ($j\in\operatorname{supp}(\bar{\vect{\theta}}^h)\Leftrightarrow|d_j|>\lambda/2$), and confirm all predictions numerically.
\end{itemize}

%================================================
% SGD as a splitting scheme
%================================================

\section{SGD as a splitting scheme}

We first consider a simple ODE where we can apply the splitting idea and the corresponding minimization problem. The best example to start from is simple ODE with right-hand-side, consisting of two summands:
\begin{equation}
    \frac{d \vect{\theta}}{d t} = - \frac{1}{2} \left( g_1(\vect{\theta}) + g_2(\vect{\theta})\right)
    \label{strang:gradientflow}
\end{equation}
Suppose, we want to find the solution $\vect{\theta}(h)$ of \eqref{strang:gradientflow} via integrating it on the small timestep $h$. The first order splitting scheme defined by solving first $\frac{d \vect{\theta}}{d t} = - \frac{1}{2} g_1(\vect{\theta}), \quad \vect{\theta}(0) = \vect{\theta}_0$ with exact solution $\vect{\theta}_1(h)$ at the moment $h$, followed by $\frac{d \vect{\theta}}{d t} = - \frac{1}{2} g_2(\vect{\theta}), \quad \vect{\theta}(0) = \vect{\theta}_1(h)$ with exact solution $\vect{\theta}_2(h)$ at the moment $h$. Thus, the first order approximation could be written as a combinations of both solutions $\vect{\theta}^I(h) = \vect{\theta}_2(h) \circ \vect{\theta}_1(h) \circ \vect{\theta}_0$.

It is interesting to study how the pure splitting scheme \cite{marchuk1968some,strang1968construction} corresponds to the SGD approach. For this purpose, we consider an illustrative example of Gradient Flow equation \ref{strang:simple_GF}, where the right-hand side of ODE is just the sum of operators acting on $\vect{\theta}$, which allows us to apply splitting scheme approximation directly.
\begin{equation}
\label{strang:simple_GF}
\frac{d \vect{\theta}}{d t} = -\frac{1}{2} \sum\limits_{i=1}^2 \nabla f_i (\vect{\theta}) = - \frac{1}{2} \nabla f_1 (\vect{\theta}) -\frac{1}{2}  \nabla f_2 (\vect{\theta})
\end{equation}
\begin{table}[h!]
\caption{The table describes the correspondence between splitting scheme for discretized Gradient Flow ODE and epoch of SGD}
\resizebox{\textwidth}{!}{
\begin{tabular}{cccc}
\toprule
\textbf{Splitting step} & \textbf{Euler discretization} & \textbf{SGD Epoch} & \textbf{First-order splitting} \\
\midrule
$\frac{d \vect{\theta}}{d t} = -\frac{1}{2}\nabla f_1(\vect{\theta})$ & $\tilde{\vect{\theta}}_{I} = \vect{\theta}_0 - \frac{h}{2}\nabla f_1 (\vect{\theta}_0) $&$\tilde{\vect{\theta}}_{SGD} = \vect{\theta}_0 - h \nabla f_1 (\vect{\theta}_0) $&$\tilde{\vect{\theta}}_{I} = \vect{\theta}_0 - \frac{h}{2}\nabla f_1 (\vect{\theta}_0)$ \\
$\frac{d \vect{\theta}}{d t} = -\frac{1}{2}\nabla f_2(\vect{\theta}) $&$\vect{\theta}_{I} = \tilde{\vect{\theta}}_{I} - \frac{h}{2}\nabla f_2 (\tilde{\vect{\theta}}_{I}) $&$\vect{\theta}_{SGD} = \tilde{\vect{\theta}}_{SGD} - h \nabla f_2 (\tilde{\vect{\theta}}_{SGD}) $&$\vect{\theta}_{I} = \tilde{\vect{\theta}}_{I} - \frac{h}{2}\nabla f_2 (\tilde{\vect{\theta}}_{I})$ \\
\bottomrule
\end{tabular}
}
\end{table}

Thus, we can conclude, that \textit{one epoch of SGD is just the splitting scheme for the discretized Gradient Flow ODE with $2 \cdot h$ step size ($m \cdot h$ in case of $m$ batches)}

Indeed, in SGD we go in the direction of the batch gradient, which stands for the Euler discretization of batch gradient flow ODE or \emph{local ODE}. This idea gives additional intuition on the method. Given information about the Euler scheme limitation (first-order accuracy, stability issues), we propose to solve each local problem more precisely.

\section{Optimization step with ODE solver}

We propose to integrate local problem more precisely instead of Euler step in SGD. Solution of the local ODE problem involves replacing gradient in the right-hand side of gradient flow ODE \ref{strang:gradientflow} with batch gradient version. In our experiments the explicit Runge-Kutta method \cite{dormand1980family,shampine1986some} was used via scipy \cite{2020SciPy-NMeth} function odeint.

\begin{table}[h!]
\caption{The table presents ODE, which we need to solve at each step of the algorithm. The last column shows the ODE, which is needed to be solved at each iteration of the algorithm for each given problem.}
\resizebox{\textwidth}{!}{
\begin{tabular}{cccc}
\toprule
\textbf{Problem} & \textbf{Loss function} & \textbf{Batch gradient} & \textbf{Initial local ODE} \\
\midrule
Linear Least Squares & $f(\vect{\theta}) = \frac{1}{n}\sum\limits_{i=1}^m\Vert X_i \vect{\theta} - \vect{y_i} \Vert_2^2$ & $\frac{1}{b}X_i^\top( X_i \vect{\theta} - \vect{y_i})$ & $\frac{d \vect{\theta}}{d t} = - \frac{1}{n} X_i^\top( X_i \vect{\theta} - \vect{y_i})$ \\
Binary logistic regression & $\begin{aligned}f(\vect{\theta}) = -\frac{1}{n} \sum_{i=1}^n\left(y_i \ln \sigma(\vect{\theta}^\top\vect{x_i})  \right.&+ \\ \left.+ (1-y_i) \ln \left(1-\sigma(\vect{\theta}^\top\vect{x_i})\right)\right)&\end{aligned}$ & $\frac{1}{b}X_i^\top\left( \sigma\left(X_i \vect{\theta}\right) - \vect{y_i}\right)$ & $\frac{d \vect{\theta}}{d t} = - \frac{1}{n} X_i^\top\left( \sigma\left(X_i \vect{\theta}\right) - \vect{y_i}\right)$ \\
One FC Layer + softmax & $f(\Theta) = 
-\frac{1}{n} \sum\limits_{i=1}^n\log\left(\frac{\vect{y_i}^\top e^{\Theta^\top \vect{x_i}}}{\vect{1}^\top e^{\Theta^\top \vect{x_i}}}\right)$ & $ \frac{1}{b} X_i^\top\left(s(\Theta^\top X_i^\top) - Y_i \right)^\top$ & $\frac{d \Theta}{d t} = - \frac{1}{n} X_i^\top\left(s(\Theta^\top X_i^\top) - Y_i \right)^\top$\\
\bottomrule
\end{tabular}
}
\end{table}

\begin{algorithm}[h!]
\SetAlgoLined
$\vect{\theta}_0$ - initial parameter; $b$ - batch size; $\alpha$ - learning rate; $m$- total number of batches

$h := \alpha m$

$t := 0$

\For{$k = 0,1, \ldots$}{
    \For{$i = 1,2, \ldots, m$}{
        Formulate local ODE problem $\mathcal{P}_i^k$

        $\vect{\theta}_{t+1} = $ integrate $\mathcal{P}_i^k$ given an initial value $\vect{\theta}(0) = \vect{\theta}_t$ to the step h

        $t := t+1$ 
        }
    }
\caption{Splitting optimization}
\end{algorithm}

Typical machine learning problems involves dealing with mini-batch of size $b$, which is often less, than the number of trainable parameters $p$, which allows us to reduce dimensionality of the dynamic system via $QR$ decomposition of each batch data matrix $X_i^\top = Q_i R_i$ (see details in the Appendix) and substitution $\vect{\eta}_i = Q_i^\top \vect{\theta}$. Note, that $QR$ decomposition is only needed to be performed once before the training.

\begin{table}[h!]
\caption{The table shows initial local ODE and paired $\mathcal{P}_i^k$. Note, that $\vect{\eta}_i \in \mathbb{R}^b$ , while $\vect{\theta} \in \mathbb{R}^p$}
\resizebox{\textwidth}{!}{
\begin{tabular}{llc}
\toprule
\multicolumn{1}{c}{\textbf{Initial local ODE}} & \multicolumn{1}{c}{$\mathcal{P}_i^k$} & \textbf{Integration} \\
\midrule
$\frac{d \vect{\theta}}{d t} = - \frac{1}{n} X_i^\top( X_i \vect{\theta} - \vect{y_i})$ & $\frac{d \vect{\eta_i}}{d t} = - \frac{1}{n} R_i\left(R_i^\top \vect{\eta_i} - \vect{y_i}\right), \vect{\eta_i} = Q_i^\top\vect{\theta}$ & analytical \\
$\frac{d \vect{\theta}}{d t} = - \frac{1}{n} X_i^\top\left( \sigma\left(X_i \vect{\theta}\right) - \vect{y_i}\right)$ & $\frac{d \vect{\eta_i}}{d t} = - \frac{1}{n} R_i\left(\sigma\left(R_i^\top \vect{\eta_i}\right) - \vect{y_i}\right), \vect{\eta_i} = Q_i^\top\vect{\theta}$ & \texttt{odeint} \\
$\frac{d \Theta}{d t} = - \frac{1}{n} X_i^\top\left(s(\Theta^\top X_i^\top) - Y_i \right)^\top$ & $\frac{d H_i}{dt} = - \frac{1}{n} R_i(s(H_i^\top R) - Y_i)^\top, H_i = Q_i^\top \Theta $& \texttt{odeint} \\
\bottomrule
\end{tabular}
}
\end{table}

There is an analytical solution for each local ODE in linear least squares case:

\begin{restatable}{theorem}{llsls}\label{strang:LLS_local_solution} For any matrix $\vect{x_i} \in \mathbb{R}^{b \times p}, b \leq p, \text{rank}X_i = b$, any vector of right-hand side $\vect{y_i} \in \mathbb{R}^{b}$ and initial vector of parameters $\vect{\theta}_0$, there is a solution of the $\frac{d \vect{\theta}}{d t} = - \frac{1}{n} X_i^\top( X_i \vect{\theta} - \vect{y_i})$, given by formula:
\begin{equation}\label{strang:LLS_local_solution_formula}
\vect{\theta}(h) = Q_i e^{-\frac{1}{n}R_iR_i^\top h} \left( Q_i^\top \vect{\theta}_0 - R_i^{-\top}\vect{y_i}\right) + Q_iR_i^{-\top}\vect{y_i} + (I - Q_iQ_i^\top)\vect{\theta}_0,
\end{equation}
where $Q_i \in \mathbb{R}^{p \times b}$ and $R_i \in \mathbb{R}^{b \times b}$ stands for the $QR$ decomposition of the matrix $\vect{X_i}^\top$, $\vect{X_i}^\top = Q_i R_i$.
\end{restatable}

It is interesting to mention, that the splitting approach immediately leads to the Kaczmarz \cite{kaczmarz1937method, strohmer2009randomized, gower2015randomized} method for solving linear system in the same setting with unit batch size.

\begin{equation}
\label{strang:splitting_limit_kaczmarz}
\lim_{h \to \infty} \vect{\theta}(h) = \frac{\left(y_i - \vect{x_i}^\top\vect{\theta}_0 \right)}{\|\vect{x_i}\|^2} \vect{x_i} + \vect{\theta}_0,
\end{equation}

which is exact formula for Kaczmarz method for solving linear system. This result correlates with the statements of \cite{needell2014stochastic}, but provides us with a new sense of similarity between SGD and Kaczmarz method.

%================================================
% Results
%================================================

\section{Results}

In this section, we describe the experimental setting. The majority of computations were performed on the NVIDIA DGX-2 cluster with 80 CPUs and 512 Gb RAM. We restricted the number of CPU usage per each experiment with an upper limit of 5 CPUs per experiment. All time measurements were done with the time library for Python. All experiments were done with the fixed random seed for reproducibility. For each experiment we performed 30 runs with random initialization and plotted trend line with the standard deviation.

\textbf{Linear Least Squares}  Both random and the real linear systems were tested. For random linear system (\texttt{random lls}) we generated $10000\times500$ matrix with additive Gaussian noise of magnitude $0.01$. Presented figures correspond to the batch size equals to $20$. The real linear system (\texttt{tom lls}) is the standard tomography data from AIRTools II \cite{hansen2018air}. Solution of the linear system is the $50 \times 50$ image reconstructed from solving $12780 \times 2500$ linear system. Presented figures correspond to the batch size equals to $60$. Relative error $10^{-3}$ was used as the stopping criterion.

\textbf{Binary Logistic Regression}  (\texttt{logreg}) In our experiments we used two classes from MNIST \cite{lecun1998gradient} dataset, which corresponds to the $0$ and $1$ digits. The size of the batch for presented figure is $50$. Test error $0.001$ was used as the stopping criterion.

\textbf{Softmax Logistic Regression} (\texttt{softmax}) We took Fashion MNIST \cite{xiao2017fashion} dataset with $60000$ grayscale pictures from $10$ classes. Each example is $28 \times 28$ image. The size of the batch for presented figure is $64$. Test error $0.25$ was used as the stopping criterion.

On the figures below we have two labels: \texttt{SGD} and \texttt{Splitting}, which stands for batch stochastic gradient descent and proposed algorithm. We use different constant learning rates to perform our experiments. All the learning rates tested for both algorithms. Lack of point of one algorithm on the graph means reaching the limit of iterations without achieving the termination rule.

\begin{figure}[h!]
    \begin{subfigure}[b]{0.25\textwidth}
            \centering
            \includegraphics[width=\linewidth]{LLS_iter.pdf}
            \caption{{\small \texttt{Random LLS}}}
    \end{subfigure}%
    \begin{subfigure}[b]{0.25\textwidth}
            \centering
            \includegraphics[width=\linewidth]{LLS_tom_iter.pdf}
            \caption{{\small \texttt{Tom LLS}}}
    \end{subfigure}%
    \begin{subfigure}[b]{0.25\textwidth}
            \centering
            \includegraphics[width=\linewidth]{Logreg_iter.pdf}
            \caption{{\small \texttt{LogReg}}}
    \end{subfigure}%
    \begin{subfigure}[b]{0.25\textwidth}
            \centering
            \includegraphics[width=\linewidth]{{Softmax_fashion_mnist_iter_err0.25}.pdf}
            \caption{{\small \texttt{Softmax}}}
    \end{subfigure}
\end{figure}

As it is expected, SGD diverges starting from some value of learning rate, which is specific for each problem. While we can see comparative robustness of the proposed splitting optimization approach.

\begin{figure}[h!]
    \begin{subfigure}[b]{0.25\textwidth}
            \centering
            \includegraphics[width=\linewidth]{LLS_time.pdf}
            \caption{{\small \texttt{Random LLS}}}
    \end{subfigure}%
    \begin{subfigure}[b]{0.25\textwidth}
            \centering
            \includegraphics[width=\linewidth]{LLS_tom_time.pdf}
            \caption{{\small \texttt{Tom LLS}}}
    \end{subfigure}%
    \begin{subfigure}[b]{0.25\textwidth}
            \centering
            \includegraphics[width=\linewidth]{Logreg_time.pdf}
            \caption{{\small \texttt{LogReg}}}
    \end{subfigure}%
    \begin{subfigure}[b]{0.25\textwidth}
            \centering
            \includegraphics[width=\linewidth]{{Softmax_fashion_mnist_time_err0.25}.pdf}
            \caption{{\small \texttt{Softmax}}}
    \end{subfigure}
\end{figure}


%================================================
% Related work
%================================================

\section{Related work}
Our work sits at the intersection of four research threads: continuous-time analysis of optimisation algorithms, operator splitting for differential equations, proximal gradient methods for sparse optimisation, and randomised linear algebra.

\textbf{Continuous-time analysis of optimisation.}
In \cite{su2014differential} authors introduced a second-order ODE that is equivalent (in the limit sense) to gradient descent with Nesterov momentum \cite{nesterov1983method}. Generalization of these ideas was presented in \cite{wibisono2016variational} with an arbitrary polynomial acceleration using the same parameter in ODE. General overview of the interplay between continuous-time and discrete-time points of view on dynamical systems and iterative optimisation methods is covered in \cite{helmke2012optimization} and \cite{evtushenko1994stable}.

\textbf{Splitting methods in numerical analysis.}
Operator splitting is a classical technique in numerical analysis for solving differential equations with composite structure. A comprehensive 161-page survey by Blanes, Casas and Murua \cite{blanes2024splitting} provides a unified treatment of splitting and composition methods with embedded error estimators, covering Lie--Trotter, Strang, and higher-order schemes. Their focus is on PDE applications (e.g.\ advection-diffusion, Hamiltonian systems) and geometric integration. Our work is complementary: we study operator splitting in the context of \emph{machine learning optimisation}, where the composite structure arises from the finite-sum objective \eqref{strang:finitesum} and the sub-problems admit closed-form solutions (LLS, softmax, LASSO). In contrast to the PDE setting, we are primarily interested in asymptotic bias and convergence properties of the resulting iterative algorithms rather than short-time approximation accuracy.

\textbf{Splitting for deep learning architectures.}
A different line of work uses splitting as an architectural principle for neural networks rather than an optimisation tool. Huang and Yin \cite{huang2022dosnet} propose Deep Operator-Splitting Networks (DOSnet), where each layer corresponds to one sub-step of a splitting scheme applied to a learned PDE or optimisation operator. Their goal is network expressiveness and interpretability, not convergence analysis of the training procedure. Our work focuses exclusively on the training dynamics of standard ML models (linear, logistic, LASSO) and makes no assumptions on network architecture.

\textbf{Proximal gradient methods and sparse optimisation.}
The LASSO problem \cite{tibshirani1996lasso} and its proximal gradient solvers are a well-studied area in optimisation. ISTA and its accelerated variant FISTA \cite{beck2009fast} are the standard first-order algorithms for composite ($\ell_2 + \ell_1$) objectives, with $O(1/K)$ and $O(1/K^2)$ convergence rates, respectively. A comprehensive treatment of proximal algorithms is given in \cite{parikh2014proximal}. The present work reinterprets ISTA as the Lie--Trotter splitting of the gradient-flow ODE for LASSO, and derives a novel result on the asymptotic bias introduced by the second-order (Strang) splitting, which has no counterpart in the proximal methods literature.

\textbf{Kaczmarz and randomised coordinate methods.}
The connection between SGD and Kaczmarz (ART) methods for linear systems has been studied in \cite{strohmer2009randomized}, which established exponential convergence of randomised Kaczmarz. We recover this algorithm as an exact limit case of the unit-batch splitting scheme for linear least-squares (see~\eqref{strang:splitting_limit_kaczmarz}), providing a new perspective via ODE splitting.

%================================================
% Conclusions
%================================================

\section{Conclusions}

We have presented a unified view of stochastic gradient methods through the lens of
operator splitting for gradient-flow ODEs.  Our main contributions are threefold.

\textbf{Splitting perspective on SGD.}
We showed that one epoch of SGD with mini-batches is precisely a first-order (Lie--Trotter)
splitting scheme applied to the full gradient-flow ODE with an effective step size
$m \cdot h$ for $m$ batches.  This connection clarifies the role of the learning rate
and batch size in the continuous-time limit, and naturally motivates replacing the
Euler-type integration of local sub-problems by higher-order ODE solvers.
Empirically, the resulting \emph{splitting optimization} scheme is considerably
more robust to the learning-rate choice than vanilla SGD across linear least squares,
logistic regression, and softmax regression benchmarks.

\textbf{Kaczmarz as a limit case.}
For the linear least-squares problem with unit-batch SGD, we established that the
splitting scheme converges to the Kaczmarz (ART) iterates in the exact-arithmetic
limit.  This sheds new light on the well-studied Kaczmarz method as the canonical
first-order splitting discretisation for linear systems.

\textbf{Splitting order matters for sparse optimisation (LASSO).}
A detailed analysis of the LASSO problem reveals a qualitative difference between
the first-order (Lie--Trotter / ISTA) and second-order (Strang-ISTA) splitting:
\begin{itemize}
  \item \textbf{LT/ISTA} converges to the exact minimiser $\vect{\theta}^*$ with zero
        asymptotic bias (Theorem~\ref{thm:lt_lasso}).
  \item \textbf{Strang-ISTA} converges to a biased fixed point $\bar{\vect{\theta}}^h$
        satisfying $F(\bar{\vect{\theta}}^h) - F^* = C_D h^2\lambda^2 + O(h^3)$,
        where $C_D \geq \|A\operatorname{sign}(\vect{\theta}^*)\|^2/8$
        (Theorem~\ref{thm:strang_ista_convergence}, Remark~\ref{rem:cd_formula}).
\end{itemize}
Hence, despite the theoretical second-order accuracy advantage of the Strang scheme
for smooth ODEs, the first-order splitting is preferable for the composite
($\ell_2 + \ell_1$) objective: the non-smooth $\ell_1$ term introduces a fixed-point
shift proportional to $h\lambda$ that the Strang scheme cannot eliminate.
Numerical experiments directly confirm the $h^2$-scaling of the bias floor and the
component-wise fixed-point shift formula~\eqref{eq:strang_shift}.
We additionally propose Conjecture~\ref{conj:strang_support}, which characterises the
support of the Strang-ISTA fixed point via the threshold rule
$j\in\operatorname{supp}(\bar{\vect{\theta}}^h)\Leftrightarrow|d_j|>\lambda/2$, and
provide numerical evidence consistent with it across a range of step sizes.

\textbf{Future directions.}
Several open questions remain.  Proving Conjecture~\ref{conj:strang_support} and
obtaining the exact formula for $C_D$ --- which depends on the support structure of
the Strang fixed point $\bar{\vect{\theta}}^h$, including off-support activations
--- are the most immediate theoretical challenges (Remark~\ref{rem:cd_formula}).
Extending the splitting analysis to general composite objectives $f + g$, combining
the splitting perspective with variance-reduction techniques (SVRG, SARAH), and
studying higher-order splitting (Yoshida--Suzuki) for structured ML problems are
natural directions.

%================================================
% References
%================================================

\newpage
\bibliographystyle{elsarticle-num}
\bibliography{biblio}

%================================================
% APPENDIX
%================================================

\appendix

%================================================
% SPLITTING GLOBAL ERROR
%================================================

\section{Upper bound on the global splitting error}

Suppose, that we have only two batches, and the problem \eqref{strang:LLS} is consistent, i.e. there exists an exact solution $\vect{\theta}_*$ such as $X \vect{\theta}_* = \vect{y}$. The GD flow has the form
\begin{equation}\label{strang:model1}
\begin{split}
    \frac{d \vect{\theta}}{d t} &= -X^{\top} (X \vect{\theta} - \vect{y}) = -X^{\top} X(\vect{\theta} - \vect{\theta}_*) =\\ &= -(X_1^{\top} X_1 + X^{\top}_2 X_2)(\vect{\theta} - \vect{\theta}_*),
\end{split}
\end{equation}
i.e. the splitting scheme corresponds to a linear operator splitting

\begin{equation*}
A = A_1 + A_2, \; A = -X^{\top} X, \; A_i = -X^{\top}_i X_i, \; i = 1, 2.
\end{equation*}

Both $A_1$ and $A_2$ are symmetric non-negative definite matrices. Without loss of generality, we can assume that $\vect{\theta}_* = 0$,

Suppose that the rank of $A$ is $r_1$ and the rank of $A_2$ is $r_2$. Then, we can write them as
\begin{equation*}
A_i = Q_i B_i Q^*_i,
\end{equation*}
where $Q_i$ is an $N \times r_i$ matrix with orthonormal columns.  The following Lemma gives the representation of the matrix exponents of such matrices.
\begin{lemma}\label{strang:lemexp}
Let $A = Q B Q^*,$ where $Q$ is an $N \times r$ matrix with orthonormal columns, and $B$ is an $r \times r$ matrix. Then,
\begin{equation}\label{strang:lrexp}
    e^{t A}  = (I - QQ^*) + Q e^{t B} Q^*.
\end{equation}
\end{lemma}
To prove \eqref{strang:lrexp} we note that
\begin{equation*}
\begin{split}
e^{t A} &= \sum_{k=0}^{\infty} \frac{t^k A^k}{k!} = \sum_{k=0}^{\infty} \frac{t^k Q B^k Q^*}{k!} = \\ &= I - QQ^* + QQ^* + Q \sum_{k=1}^{\infty} \frac{t^k B^k}{k!} Q^* = \\ &= (I - QQ^*) + Q e^{t B} Q^*.
\end{split}
\end{equation*}

\begin{lemma}
    \label{strang:lemupper_2}
    Let $A_1, A_2 \in \mathbb{S}^p_{+}$ be the square negative semidefinite matrices, that don't have full rank, i.e. $\operatorname{rank}{A_1} \leq p$ and $\operatorname{rank}{A_2} \leq p$. While the sum of those matrices has full rank, i.e. $A = A_1 + A_2, \operatorname{rank}{A} = p$. Then, the global upper bound error will be written as follows:

    \begin{equation}\label{strang:lemupper}
        \lim_{t \to \infty}\| e^{A_2t}e^{A_1t} - e^{At}\| = \|(I - Q_2Q_2^*)(I - Q_1Q_1^*)\|
    \end{equation}
\end{lemma}
\begin{proof}
    The proof is straightforward. We will use the low rank matrix exponential decomposition from the Lemma \ref{strang:lemexp}
    $$
    e^{A_it} = \Pi_i + Q_i e^{B_it} Q_i^*, \text{where } \Pi_i = I - Q_iQ_i^*; i = 1,2
    $$
    \begin{align*}
    &\lim_{t \to \infty}\| e^{A_2t}e^{A_1t} - e^{At}\| = \\ 
    &= \lim_{t \to \infty}\| (\Pi_2 + Q_2 e^{B_2t} Q_2^*)(\Pi_1 + Q_1 e^{B_1t} Q_1^*) - e^{At}\| = \\
    &= \lim_{t \to \infty}\| \Pi_2\Pi_1 + Q_1 e^{B_1t} Q_1^*\Pi_2 + \Pi_1Q_2 e^{B_2t} Q_2^* + \\
    &+  Q_1 e^{B_1t} Q_1^* Q_2 e^{B_2t} Q_2^* - e^{At}\| =\\&= \Pi_2 \Pi_1
    \end{align*}
    Since all matrices $B_1, B_2, A$ are negative all the matrix exponentials are decaying: $\|e^{At}\|\leq e^{t\mu (A)}\, \forall t\geq 0$, where $\mu(A) = \lambda_{max} \left( \frac{A + A^\top}{2}\right)$ - the logarithmic norm.
\end{proof}

\begin{figure}[h!]
    \begin{subfigure}[t]{0.49\textwidth}
            \centering
            \includegraphics[width=\linewidth]{upper_bound_2.pdf}
            \caption{{\small Global error of the splitting scheme. Initial random full rank matrix $X \in \mathbb{R}^{100 \times 100}$ was splitted by rows. $X_1, X_2 \in \mathbb{R}^{50 \times 100}$. Target matrices were obtained the following way: $A_1 = -X_1^*X_1, A_2 = -X_2^*X_2, A = -X^*X$. So $A_1, A_2$ are negative and lacking full rank, while $A = A_1 + A_2$ has full rank.}}
            \label{strang:fig:upper_bound_2}
    \end{subfigure}%
    \hfill
    \begin{subfigure}[t]{0.49\textwidth}
            \centering
            \includegraphics[width=\linewidth]{upper_bound_many.pdf}
            \caption{{\small Global upper bound on the splitting scheme in case of $40$ summands in the right-hand side.}}
            \label{strang:fig:upper_bound_many}
    \end{subfigure}
\end{figure}

The graph presented on the Figure \ref{strang:fig:upper_bound_2} describes the global splitting error and the derived upper bound compared to the existing bound from \cite{sheng1994global} for the case of $2$ summands in the right-hand side. One can easily see significant difference between existing global upper bounds for that case and the derived upper bound.

\begin{theorem}\label{strang:theorem_uppbound}
    Let $A_1, A_2, \ldots, A_b \in \mathbb{S}^p_{+}$ be the square negative semidefinite matrices, that don't have full rank, i.e. $\operatorname{rank}{A_i} \leq p, \;\forall i = 1, \ldots, b$. While the sum of those matrices has full rank, i.e. $A = \sum\limits_{i=1}^b A_i, \operatorname{rank}{A} = p$. Then, the global upper bound error will be written as follows:

    \begin{equation}\label{strang:global_error_upper_bound}
        \lim_{t \to \infty}\| e^{A_bt} \cdot \ldots \cdot e^{A_1t} - e^{At}\| = \left\|\prod\limits_{i=1}^b \Pi_{b-i+1}\right\|,
    \end{equation}
    where $\Pi_i = I - Q_iQ_i^*$ and $A_i = Q_iB_iQ_i^*$ and $Q_i$ is a matrix with orthonormal columns. 
\end{theorem}

The graph on the Figure \ref{strang:fig:upper_bound_many} shows empirical validity of the presented upper bound.


%================================================
% Proofs
%================================================

\section{Proofs}

\llsls*

\begin{proof}
Given $X_i^\top = Q_i R_i$, we have $(I - Q_iQ_i^\top)X_i^\top = 0$. Note, that $Q_i$ is left unitary matrix, i.e. $Q_i^\top Q_i = I$. 
\begin{align}
\nonumber \frac{d \vect{\theta}}{d t} &= - \frac{1}{n}X_i^\top(X_i \vect{\theta} - \vect{y_i}) \\
\nonumber (I - Q_iQ_i^\top)\frac{d \vect{\theta}}{d t} &= 0 \\
\nonumber \frac{d \vect{\theta}}{d t} &= Q_i\frac{d (Q_i^\top\vect{\theta})}{d t} \quad Q_i^\top \vect{\theta} = \vect{\eta_i}\\ 
\nonumber \frac{d \vect{\theta}}{d t} &= Q_i\frac{d \vect{\eta_i}}{d t} \quad \text{integrate from $0$ to $h$}\\ 
\label{strang:lls_theorem_theta_from_eta}\vect{\theta}(h) &= Q_i \left(\vect{\eta_i}(h) - \vect{\eta_i}(0) \right) + \vect{\theta}_0
\end{align}

On the other hand:
\begin{align}\nonumber
\frac{d \vect{\eta_i}}{d t} &= Q_i^\top\frac{d \vect{\theta}}{d t} =  - \frac{1}{n} Q_i^\top X_i^\top(X_i \vect{\theta} - \vect{y_i}) = \\ 
\nonumber&= - \frac{1}{n} Q_i^\top  Q_i R_i( R_i^\top Q_i^\top \vect{\theta} - \vect{y_i}) =\\
&= - \frac{1}{n} \left( R_iR_i^\top \vect{\eta_i} - R_i \vect{y_i}\right) \label{strang:lls_theorem_eta_from_theta}
\end{align}

Consider the moment of time $t = \infty$. $\frac{d \vect{\eta_i}}{d t} = 0$, since $\exists \vect{\theta}^*, Q_i^\top \vect{\theta}^* = \vect{\eta_i}^*$. Also consider \eqref{strang:lls_theorem_eta_from_theta}:

\begin{equation}\label{strang:lls_theorem_eta_star}
\begin{split}
\frac{d \vect{\eta_i}}{d t} = 0 &= - \frac{1}{n} \left( R_iR_i^\top \vect{\eta_i}^* - R_i \vect{y_i}\right) \\ R_i \vect{y_i} &= R_iR_i^\top \vect{\eta_i}^*
\end{split}
\end{equation}

Now we look at the \eqref{strang:lls_theorem_eta_from_theta} with the replacement, given in \eqref{strang:lls_theorem_eta_star}:

\begin{align*}
\frac{d \vect{\eta_i}}{d t} &= - \frac{1}{n} \left( R_iR_i^\top \vect{\eta_i}- R_iR_i^\top \vect{\eta_i}^*\right) \\
\frac{d \vect{\eta_i}}{d t} &= - \frac{1}{n} R_iR_i^\top \left( \vect{\eta_i}- \vect{\eta_i}^*\right) \qquad \text{integrate from $0$ to $h$} \\
\vect{\eta_i}(h) - \vect{\eta_i}^* &= e^{- \frac{1}{n} R_iR_i^\top h} (\vect{\eta_i}(0) - \vect{\eta_i}^*) \\ 
&\text{while } \vect{\eta_i}^* = R_i^{-\top}\vect{y_i}, \vect{\eta_i}(0) = Q_i^\top \vect{\theta}_0 \\
\vect{\eta_i}(h) &= e^{- \frac{1}{n} R_iR_i^\top h} (Q_i^\top \vect{\theta}_0 - R_i^{-\top}\vect{y_i}) + R_i^{-\top}\vect{y_i} 
\end{align*}

Using \eqref{strang:lls_theorem_theta_from_eta} we obtain the target formula
\begin{equation*}
\begin{split}
\vect{\theta}(h) &= Q_i e^{-\frac{1}{n}R_iR_i^\top h} \left( Q_i^\top \vect{\theta}_0 - R_i^{-\top}\vect{y_i}\right) + \\ &+ Q_iR_i^{-\top}\vect{y_i} + (I - Q_iQ_i^\top)\vect{\theta}_0,
\end{split}
\end{equation*}

\end{proof}

%================================================
% Applications
%================================================

\section{Applications}
\subsection{Linear least squares}
\subsubsection{Problem}
Let $f_i(\vect{\theta}) = \Vert \vect{x_i}^{\top} \vect{\theta} - y_i \Vert^2$, then problem \eqref{strang:finitesum} is the linear least squares problem, which can be written as
\begin{equation}\label{strang:LLS}
   f(\vect{\theta}) = \frac{1}{n}\Vert X \vect{\theta} - \vect{y} \Vert_2^2  = \frac{1}{n}\sum\limits_{i=1}^s\Vert X_i \vect{\theta} - \vect{y_i} \Vert_2^2\to \min_{\vect{\theta} \in \mathbb{R}^p},
\end{equation}
where $X \in \mathbb{R}^{n \times p}$ and $\vect{y} \in \mathbb{R}^p$ and the second part of the equation stands for $s$ mini-batches with size $b$ regrouping ($b \cdot s = n$): $X_i \in \mathbb{R}^{b \times p}, \vect{y_i} \in \mathbb{R}^{b}$

\begin{equation}\label{strang:LLS_grad}
\nabla_\theta f(\vect{\theta}) = \nabla f(\vect{\theta}) = \frac{1}{n}\sum\limits_{i=1}^s X_i^\top(X_i \vect{\theta} - \vect{y_i})
\end{equation}

The gradient flow equation will be written as follows:
\begin{equation}\label{strang:LLS_GF}
\frac{d \vect{\theta}}{d t} = - \frac{1}{n}\sum\limits_{i=1}^s X_i^\top( X_i \vect{\theta} - \vect{y_i})
\end{equation}

\subsubsection{Exact solution of the local problem}
Theorem \ref{strang:LLS_local_solution} gives us explicit formula for the local solution:
$$
\vect{\theta}(h) = Q_i e^{-\frac{1}{n}R_iR_i^\top h} \left( Q_i^\top \vect{\theta}_0 - R_i^{-\top}\vect{y_i}\right) + Q_iR_i^{-\top}\vect{y_i} + (I - Q_iQ_i^\top)\vect{\theta}_0
$$

\subsubsection{Kaczmarz as the limit case of splitting}
Kaczmarz method \cite{kaczmarz1937method}, \cite{strohmer2009randomized}, \cite{gower2015randomized} is a well-known iterative algorithm for solving linear systems
It is interesting to mention, that splitting approach immediately leads to the Kaczmarz method for solving linear system in the same setting with unit batch size.

When the batch size is equal to one, we need to do $n$ QR decompositions for each transposed batch matrix, which is just column vector $\vect{x_i}$ in our case:

\begin{equation}
\vect{x_i} = \vect{q_i} \vect{r_i} = \underset{\vect{q_i} }{\frac{\vect{x_i}}{\|\vect{x_i}\|}} \underset{\vect{r_i}}{\vphantom{\frac{\vect{x_i}}{\|\vect{x_i}\|}} \|\vect{x_i}\|}
\end{equation}

Now, we need to use \eqref{strang:LLS_local_solution_formula} to derive analytic local solution in that case:

\begin{equation*}
\begin{split}
\vect{\theta}(h) &= \frac{\vect{x_i}}{\|\vect{x_i}\|} e^{-\frac{\|\vect{x_i}\|^2 h}{n}} \left( \frac{\vect{x_i}^\top}{\|\vect{x_i}\|} \vect{\theta}_0 - \frac{y_i}{\|\vect{x_i}\|}\right) + \\ &+ \frac{\vect{x_i}}{\|\vect{x_i}\|^2}y_i + \left(I - \frac{\vect{x_i}\vect{x_i}^\top}{\|\vect{x_i}\|^2}\right)\vect{\theta}_0 = \\
&= \frac{\left(y_i -\vect{x_i}^\top\vect{\theta}_0 \right)}{\|\vect{x_i}\|^2}  \left(1 - e^{-\frac{\|\vect{x_i}\|^2 h}{n}}\right)\vect{x_i} + \vect{\theta}_0
\end{split}
\end{equation*}

It can be easily seen, that:

\begin{equation}
\label{strang:kaczmarz_derivation}
\lim_{h \to \infty} \vect{\theta}(h) = \frac{\left(y_i - \vect{x_i}^\top\vect{\theta}_0 \right)}{\|\vect{x_i}\|^2} \vect{x_i} + \vect{\theta}_0,
\end{equation}
which is exact formula for Kaczmarz method for solving linear system. This result correlates with the statements of \cite{needell2014stochastic}, but provides us with a new sense of similarity between SGD and Kaczmarz method.

\subsection{Binary logistic regression}
\subsubsection{Problem}
In this classification task then problem \eqref{strang:finitesum} takes the following form:
\begin{equation}\label{strang:LogReg}
-\frac{1}{n} \sum_{i=1}^n\left(y_i \ln \sigma(\vect{\theta}^\top\vect{x_i})  + (1-y_i) \ln (1-\sigma(\vect{\theta}^\top\vect{x_i}))\right) \to \min_{\vect{\theta} \in \mathbb{R}^p},
\end{equation}
where $\sigma(x) = \frac{1}{1 + e^{-x}}$ is the sigmoid function, while $ y_i \in \{0,1\}$ stands for the label of the object class.

\begin{equation}\label{strang:LogReg_grad}
\nabla_\theta f(\vect{\theta}) = \nabla f(\vect{\theta}) = \frac{1}{n}\sum\limits_{i=1}^n \vect{x_i}(\sigma(\vect{\theta}^\top\vect{x_i}) - y_i)
\end{equation}

The gradient flow equation will be written as follows:
\begin{equation}\label{strang:LogReg_GF}
\frac{d \vect{\theta}}{d t} = - \frac{1}{n}\sum\limits_{i=1}^n \vect{x_i}(\sigma(\vect{\theta}^\top\vect{x_i}) - y_i)
\end{equation}

Our particular interest lies in mini-batch reformulation of the given problem. We consider $s$ mini-batches with size $b$ regrouping ($b \cdot s = n$): $X_i \in \mathbb{R}^{b \times p}, \vect{y_i} \in \mathbb{R}^{b}$ and  $\sigma(\vect{x})$ stands for the element-wise sigmoid function.

\begin{equation}\label{strang:LogReg_GF_batch}
\frac{d \vect{\theta}}{d t} = - \frac{1}{n}\sum\limits_{i=1}^s X_i^\top\left(\sigma\left(X_i\vect{\theta}\right) - \vect{y_i}\right)
\end{equation}

\subsubsection{Splitting scheme and local problem}

Since we are applying splitting scheme to find the approximate solution of the \eqref{strang:LogReg_GF_batch}, each local problem should be written as follows:

\begin{equation}\label{strang:LogReg_GF_local}
\frac{d \vect{\theta}}{d t} =  - \frac{1}{n} X_i^\top\left(\sigma\left(X_i\vect{\theta}\right) - \vect{y_i}\right)
\end{equation}

Note, that this is not linear equation and cannot be solved as easy as in Theorem \ref{strang:LLS_local_solution}. However, we can apply the same technique to reduce the dimension of ODE, which is needed to be solved numerically. 

Suppose, we have $QR$ decomposition of each batch data matrix $X_i^\top = Q_i R_i$, then we can multiply both sides of \eqref{strang:LogReg_GF_local} on the $(I - Q_iQ_i^\top)$ on the left. 

\begin{align}
\nonumber (I - Q_iQ_i^\top)\frac{d \vect{\theta}}{d t} &= (I - Q_iQ_i^\top) \frac{1}{n} X_i^\top(\vect{y_i} -\sigma\left(X_i\vect{\theta}\right)) \\
\nonumber \frac{d \vect{\theta}}{d t} &= Q_i\frac{d (Q_i^\top\vect{\theta})}{d t} \quad Q_i^\top \vect{\theta} = \vect{\eta_i}\\ 
\nonumber \frac{d \vect{\theta}}{d t} &= Q_i\frac{d \vect{\eta_i}}{d t} \quad \text{integrate from $0$ to $h$}\\ 
\label{strang:logreg_theta_from_eta}\vect{\theta}(h) &= Q_i \left(\vect{\eta_i}(h) - \vect{\eta_i}(0) \right) + \vect{\theta}_0
\end{align}

On the other hand:
\begin{align}\nonumber
\frac{d \vect{\eta_i}}{d t} &= Q_i^\top\frac{d \vect{\theta}}{d t} =  - \frac{1}{n} Q_i^\top  X_i^\top(\sigma\left(X_i\vect{\theta}\right) - \vect{y_i}) = \\ 
\nonumber&= - \frac{1}{n} Q_i^\top  Q_i R_i(\sigma\left(X_i\vect{\theta}\right) - \vect{y_i}) =\\
\nonumber&= - \frac{1}{n} R_i(\sigma\left(X_i\vect{\theta}\right) - \vect{y_i})
\end{align}

Recall, that each hypothesis function depends on linear function $\vect{x_i}^\top \vect{\theta}$, which means, that in batch reformulation it is just entries of the vector $X_i \vect{\theta}$. Since we have $QR$ decomposition of $X_i^\top$, we can write: $X_i \vect{\theta} = R_i^\top Q_i^\top \vect{\theta} = R_i^\top \vect{\eta_i}$. In other words:

\begin{equation}
\label{strang:logreg_eta_ode}
\frac{d \vect{\eta_i}}{d t} = - \frac{1}{n} R_i\left(\sigma\left(R_i^\top \vect{\eta_i}\right) - \vect{y_i}\right),
\end{equation}

To sum it up, we need to solve \eqref{strang:logreg_eta_ode} (which is much simpler, than original differential equation \eqref{strang:LogReg_GF_local}), than substitute it to the
\eqref{strang:logreg_theta_from_eta} with $\vect{\eta_i}(0) = Q_i^\top \vect{\theta}_0$. Note, that matrices $Q_i$ and $R_i$ can be computed only once before the training.

\subsection{Softmax Regression}
\subsubsection{Problem}
In this classification task then problem \eqref{strang:finitesum} takes the following form:
\begin{equation}\label{strang:Softmax}
-\frac{1}{n} \sum_{i=1}^n\log\left(\frac{\vect{y_i}^\top e^{\Theta^\top \vect{x_i}}}{\vect{1}^\top e^{\Theta^\top \vect{x_i}}}\right) \to \min_{\Theta \in \mathbb{R}^{p \times K}},
\end{equation}
where $e^{\vect{x}}$ is element-wise exponential function, while $ \vect{y_i} \in \mathbb{R}^K$ stands for the one-hot encoding of the $i$-th object label.

\begin{equation}
\nabla_\Theta f(\Theta) = -\frac{1}{n} \sum_{i=1}^n\vect{x_i}\left(\vect{y_i} - \frac{ e^{\Theta^\top \vect{x_i}}}{\vect{1}^\top e^{\Theta^\top \vect{x_i}}}\right)^\top
\end{equation}

\begin{equation}
\nabla_\Theta f(\Theta) = -\frac{1}{n} \sum_{i=1}^n\vect{x_i}\left(\vect{y_i} - s\left(\Theta^\top \vect{x_i}\right)\right)^\top
\end{equation}

% \begin{equation}
% \nabla_{\theta_k} f = -\frac{1}{n} \sum_{i=1}^n\left(y_i^k - \frac{ e^{\theta_k^\top \vect{x_i}}}{\vect{1}^\top e^{\Theta \vect{x_i}}}\right)\vect{x_i}^\top
% \end{equation}

Here we use $s(\vect{x})$ as a softmax function of a vector $\vect{x}$, i.e. $s(\vect{x}) = \frac{e^{\vect{x}}}{\vect{1}^\top e^{\vect{x}}}$ .While mini-batch reformulation will take the following form:

\begin{equation}
\nabla_\Theta f(\Theta) = -\frac{1}{n} \sum_{i=1}^s X_i^\top\left(Y_i - s(\Theta^\top X_i^\top)\right)^\top,
\end{equation}

where $s(X) = \left[\begin{array}{cccc}| & | & | & | \\
s(\vect{x}_{(1)}) & s(\vect{x}_{(2)}) & \cdots & s(\vect{x}_{(b)}) \\
| & | & | & |
\end{array}\right]$ is a column-wise softmax function. Indeed, in a very similar manner to the binary logistic regression we can write down gradientflow ODE for softmax regression in a mini-batch form:

\begin{equation}
\frac{d \Theta}{d t} = - \frac{1}{n} \sum_{i=1}^s X_i^\top\left(s(\Theta^\top X_i^\top) - Y_i \right)^\top
\end{equation}

Splitting method requires the local problem, which is focused on a single minibatch:

\begin{equation}
\frac{d \Theta}{d t} = - \frac{1}{n} X_i^\top\left(s(\Theta^\top X_i^\top) - Y_i \right)^\top
\end{equation}

\begin{align}
\nonumber (I - Q_iQ_i^\top)\frac{d \Theta}{d t} &= (I - Q_iQ_i^\top) \frac{1}{n} X_i^\top(Y_i -s(\Theta^\top X_i^\top))^\top \\
\nonumber \frac{d \Theta}{d t} &= Q_i\frac{d (Q_i^\top\Theta)}{d t} \quad Q_i^\top \Theta = H_i \\ 
\nonumber \frac{d \Theta}{d t} &= Q_i\frac{d H_i }{d t} \quad \text{integrate from $0$ to $h$}\\ 
\label{strang:softmax_theta_from_eta}\Theta(h) &= Q_i \left(H_i (h) - H_i (0) \right) + \Theta_0
\end{align}

On the other hand:
\begin{align}\nonumber
\frac{d H_i }{d t} &= Q_i^\top\frac{d \Theta}{d t} =  - \frac{1}{n} Q_i^\top  X_i^\top(s(\Theta^\top X_i^\top) - Y_i)^\top = \\ 
\nonumber&= - \frac{1}{n} Q_i^\top  Q_i R_i(s(\Theta^\top X_i^\top) - Y_i)^\top =\\
\nonumber&= - \frac{1}{n} R_i(s(\Theta^\top X_i^\top) - Y_i)^\top =\\
\nonumber&= - \frac{1}{n} R_i(s(H_i^\top R) - Y_i)^\top 
\end{align}

Now we need to solve ODE of variable of the size $b \times k$, rather, than $p \times k$.

%================================================
% LASSO with Splitting Schemes
% Status: DRAFT — Theorems A and D are sketches.
%         Proofs to be finalized with I. Oseledets (April 2026).
%================================================

\section{LASSO with Splitting Schemes}
\label{sec:lasso}

We now consider the LASSO problem \cite{tibshirani1996lasso} as a concrete instance where the choice of splitting
order has a measurable impact on convergence. Let
\begin{equation}
  F(\vect{\theta}) = \frac{1}{2}\|A\vect{\theta} - \vect{y}\|^2
    + \lambda \|\vect{\theta}\|_1 \;\to\; \min_{\vect{\theta} \in \mathbb{R}^d},
  \label{eq:lasso}
\end{equation}
where $A \in \mathbb{R}^{n \times d}$, $\vect{y} \in \mathbb{R}^n$, and $\lambda > 0$.
Let $L = \|A^\top A\|$ (spectral norm), $\vect{\theta}^*$ a minimizer of~\eqref{eq:lasso},
and $0 < h \leq 1/L$.

\subsection{Lie–Trotter (Proximal Gradient / ISTA)}

The first-order (Lions–Teboulle / LT) splitting applied to the gradient-flow ODE
for~\eqref{eq:lasso} yields the standard proximal gradient step:
\begin{equation}
  \vect{\theta}_{k+1} = \mathrm{prox}_{h\lambda\|\cdot\|_1}\!\bigl(\vect{\theta}_k - h A^\top(A\vect{\theta}_k - \vect{y})\bigr),
  \label{eq:lt_step}
\end{equation}
which coincides with ISTA \cite{beck2009fast}. The following result characterises its
asymptotic bias.

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

\begin{proof}
Write $F = f + g$ where $f(\vect{\theta}) = \tfrac{1}{2}\|A\vect{\theta} - \vect{y}\|^2$
(smooth, $L$-smooth with $L = \|A^\top A\|$) and $g(\vect{\theta}) = \lambda\|\vect{\theta}\|_1$
(convex, non-smooth).

\medskip\noindent\textbf{Step 1 (Descent lemma).}
Since $f$ is $L$-smooth and convex,
\begin{equation}\label{eq:lasso_descent}
  f(\vect{u}) \;\leq\; f(\vect{\theta}_k)
    + \langle \nabla f(\vect{\theta}_k),\, \vect{u} - \vect{\theta}_k\rangle
    + \tfrac{L}{2}\|\vect{u} - \vect{\theta}_k\|^2
  \qquad \forall\,\vect{u}.
\end{equation}

\medskip\noindent\textbf{Step 2 (Subgradient at prox-step).}
The LT iterate $\vect{\theta}_{k+1} = \operatorname{prox}_{h\lambda\|\cdot\|_1}
(\vect{\theta}_k - h\nabla f(\vect{\theta}_k))$
minimises the composite quadratic
$Q_k(\vect{\theta}) := \langle\nabla f(\vect{\theta}_k), \vect{\theta}\rangle
  + g(\vect{\theta}) + \tfrac{1}{2h}\|\vect{\theta} - \vect{\theta}_k\|^2$.
By first-order optimality there exists $\vect{p}_{k+1} \in \partial g(\vect{\theta}_{k+1})$ with
\begin{equation}\label{eq:lasso_opt_cond}
  \vect{p}_{k+1} + \nabla f(\vect{\theta}_k)
    + \tfrac{1}{h}(\vect{\theta}_{k+1} - \vect{\theta}_k) = \vect{0}.
\end{equation}

\medskip\noindent\textbf{Step 3 (One-step bound).}
Apply~\eqref{eq:lasso_descent} with $\vect{u} = \vect{\theta}_{k+1}$, add $g(\vect{\theta}_{k+1})$,
and use the subgradient inequality $g(\vect{\theta}_{k+1}) \leq g(\vect{\theta}^*)
  + \langle\vect{p}_{k+1},\, \vect{\theta}_{k+1} - \vect{\theta}^*\rangle$.
Substituting $\vect{p}_{k+1}$ from~\eqref{eq:lasso_opt_cond} and applying the identity
$\langle\vect{a} - \vect{b},\, \vect{a} - \vect{c}\rangle
  = \tfrac{1}{2}(\|\vect{a} - \vect{c}\|^2 - \|\vect{b} - \vect{c}\|^2 + \|\vect{a} - \vect{b}\|^2)$
together with convexity of $f$ gives, for $h \leq 1/L$:
\begin{equation}\label{eq:lasso_one_step}
  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).
\end{equation}

\medskip\noindent\textbf{Step 4 (Telescoping).}
Summing~\eqref{eq:lasso_one_step} over $k = 0,\ldots,K-1$ and applying convexity of $F$:
\[
  F(\bar{\vect{\theta}}_K) \;\leq\; \tfrac{1}{K}\sum_{k=0}^{K-1} F(\vect{\theta}_{k+1})
  \;\leq\; F(\vect{\theta}^*) + \frac{\|\vect{\theta}_0 - \vect{\theta}^*\|^2}{2hK},
\]
which is the stated bound. Since the right-hand side tends to zero,
the LT iterate has zero algorithmic bias.
\end{proof}

\subsection{Strang Splitting (Strang-ISTA)}

The symmetric (Strang) splitting applies the proximal operator at half the step on both
sides of the gradient step:
\begin{equation}
  \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).
  \label{eq:strang_step}
\end{equation}

Unlike the LT scheme, the Strang splitting introduces an \emph{algorithmic bias}: the
iterates converge not to $\vect{\theta}^*$ but to a shifted fixed point
$\bar{\vect{\theta}}^h \neq \vect{\theta}^*$, with a bias floor proportional to $h^2\lambda^2$.

\begin{theorem}[Strang-ISTA convergence with bias floor]\label{thm:strang_ista_convergence}
  Let $\bar{\vect{\theta}}_K = \frac{1}{K}\sum_{k=0}^{K-1}\vect{\theta}_{k+1}$ be the averaged
  iterates of~\eqref{eq:strang_step}. Then
  \begin{equation}
    F(\bar{\vect{\theta}}_K) - F(\vect{\theta}^*)
    \;\leq\;
    \frac{\|\vect{\theta}_0 - \vect{\theta}^*\|^2}{2hK} + C_D \cdot h^2 \cdot \lambda^2,
    \label{eq:strang_rate}
  \end{equation}
  where $C_D \geq \|A\operatorname{sign}(\vect{\theta}^*)\|^2/8$ is a constant depending on $A$
  and $\vect{\theta}^*$ (see Remark~\ref{rem:cd_formula}).
  In contrast to Theorem~\ref{thm:lt_lasso}, the second term is an irreducible bias
  that does not vanish as $K \to \infty$.
\end{theorem}

\begin{proof}[Proof sketch]
We outline the argument; a complete proof will appear in the journal version
(planned: discussion with I.~Oseledets, April 2026).

\medskip\noindent\textbf{Step 1 (Modified objective and biased fixed point).}
Define $\tilde{f}_h(\vect{\theta}) := f(\operatorname{prox}_{h\lambda/2}(\vect{\theta}))
  + \tfrac{1}{2h}\|\operatorname{prox}_{h\lambda/2}(\vect{\theta}) - \vect{\theta}\|^2$
and $\tilde{F}_h := \tilde{f}_h + (h\lambda/2)\|\cdot\|_1$.
One verifies that the Strang-ISTA step~\eqref{eq:strang_step} is a proximal gradient
step on $\tilde{F}_h$, so the sequence $\{\vect{\theta}_k\}$ converges to the fixed point
$\bar{\vect{\theta}}^h$ satisfying $\vect{0} \in \partial \tilde{F}_h(\bar{\vect{\theta}}^h)$.
In general, $\bar{\vect{\theta}}^h \neq \vect{\theta}^*$.

\medskip\noindent\textbf{Step 2 (ISTA-type bound relative to fixed point).}
Applying the same argument as the proof of Theorem~\ref{thm:lt_lasso}
to the sequence $\{\vect{\theta}_k\}$ with $\vect{\theta}^*$ replaced by $\bar{\vect{\theta}}^h$:
\begin{equation}\label{eq:strang_to_fp}
  \tilde{F}_h(\bar{\vect{\theta}}_K) - \tilde{F}_h(\bar{\vect{\theta}}^h)
  \;\leq\; \frac{\|\vect{\theta}_0 - \bar{\vect{\theta}}^h\|^2}{2hK}.
\end{equation}

\medskip\noindent\textbf{Step 3 (Fixed-point shift and bias floor).}
Using $\tilde{F}_h(\vect{u}) = F(\vect{u}) + O(h^2\lambda^2)$ and the fixed-point shift
$\bar{\vect{\theta}}^h - \vect{\theta}^* \approx (h\lambda/2)\operatorname{sign}(\vect{\theta}^*)$
on the support of $\vect{\theta}^*$
(see Remark~\ref{rem:cd_formula}), one obtains
$F(\bar{\vect{\theta}}^h) - F(\vect{\theta}^*) \geq \|A\operatorname{sign}(\vect{\theta}^*)\|^2 h^2\lambda^2/8 + O(h^3)$;
off-support components activated at the fixed point contribute additional positive terms,
so $C_D \geq \|A\operatorname{sign}(\vect{\theta}^*)\|^2/8$.

\medskip\noindent\textbf{Step 4 (Combining).}
Substituting into~\eqref{eq:strang_to_fp} and using
$\|\vect{\theta}_0 - \bar{\vect{\theta}}^h\| \leq \|\vect{\theta}_0 - \vect{\theta}^*\| + O(h\lambda)$
yields the stated bound after absorbing lower-order terms into $C_D h^2\lambda^2$.
\end{proof}

\begin{remark}[Constant $C_D$ and Strang fixed-point shift]
\label{rem:cd_formula}
The constant $C_D$ in Theorem~\ref{thm:strang_ista_convergence} characterises the
algorithmic bias introduced by Strang splitting applied to the LASSO problem.
Numerical estimation
(for $A \in \mathbb{R}^{200 \times 20}$, $\lambda = 0.05$, varying $h \in [0.1/L, 0.5/L]$)
yields
\begin{equation}
  C_D \approx 0.96.
  \label{eq:cd_empirical}
\end{equation}

\noindent \textbf{Analytical lower bound.}
On the \emph{support} of $\vect{\theta}^*$, the fixed-point shift equals
$\bar{\vect{\theta}}^h_i - \vect{\theta}^*_i = (h\lambda/2)\operatorname{sign}(\vect{\theta}^*_i)$
exactly (for $h$ small enough that the support is preserved).
Expanding $F$ around $\vect{\theta}^*$ using this on-support shift gives
\begin{equation}
  C_D \;\geq\; \frac{\|A \cdot \operatorname{sign}(\vect{\theta}^*)\|^2}{8},
  \label{eq:cd_formula}
\end{equation}
where $\vect{\theta}^*$ is the optimal solution.
For the numerical setting above: $\|A\,\vect{s}^*\|^2 \approx 6.15$,
giving a lower bound $\|A\,\vect{s}^*\|^2/8 \approx 0.769$.

However, the Strang fixed point can activate \emph{off-support} components
(indices $i \notin \operatorname{supp}(\vect{\theta}^*)$) where the dual variable satisfies
$|[A^\top(A\bar{\vect{\theta}}^h - \vect{y})]_i| \approx \lambda/2 \cdot (1 - \epsilon)$
for some $0 < \epsilon < 1$.
These off-support activations contribute positive terms to $F(\bar{\vect{\theta}}^h) - F(\vect{\theta}^*)$
beyond the on-support lower bound~\eqref{eq:cd_formula}, so in general $C_D > \|A\,\vect{s}^*\|^2/8$.
Numerical experiments confirm $C_D \approx 0.96 > 0.769$ for the test problem
(see~\eqref{eq:cd_empirical}).

\noindent \textbf{Fixed-point shift.}
The Strang fixed point $\bar{\vect{\theta}}^h$ is shifted from $\vect{\theta}^*$ by
\begin{equation}
  \bar{\vect{\theta}}^h_i - \vect{\theta}^*_i \approx \frac{h\lambda}{2}\operatorname{sign}(\vect{\theta}^*_i),
  \qquad i \in \operatorname{supp}(\vect{\theta}^*),
  \label{eq:strang_shift}
\end{equation}
on the support of $\vect{\theta}^*$; off-support components may be non-zero.
This bias is proportional to the step size $h$ and regularisation strength $\lambda$,
and vanishes in the limit $h \to 0$ (consistent with the zero bias of LT).
The excess objective value satisfies
\begin{equation}
  F(\bar{\vect{\theta}}^h) - F(\vect{\theta}^*)
  \;\geq\; \frac{h^2\lambda^2}{8}(\vect{s}^*)^\top H\,\vect{s}^* + O(h^3),
\end{equation}
where $H = A^\top A$ and $\vect{s}^* = \operatorname{sign}(\vect{\theta}^*)$.
The exact constant $C_D$ depends on the full support structure of the Strang fixed point
$\bar{\vect{\theta}}^h$ and requires further analysis
(planned: discuss with I.~Oseledets, April 2026).
\end{remark}

\begin{conjecture}[Strang support characterisation]\label{conj:strang_support}
Let $\vect{d} = A^\top(A\vect{\theta}^* - \vect{y})$ be the dual variable at the LT
fixed point $\vect{\theta}^*$.  For sufficiently small $h > 0$, the support of the
Strang-ISTA fixed point $\bar{\vect{\theta}}^h$ satisfies
\begin{equation}\label{eq:strang_support_char}
  j \;\in\; \operatorname{supp}(\bar{\vect{\theta}}^h)
  \;\Longleftrightarrow\;
  |d_j| \;>\; \tfrac{\lambda}{2}.
\end{equation}
In particular:
\begin{itemize}
  \item all support indices $j \in \operatorname{supp}(\vect{\theta}^*)$ satisfy $|d_j| = \lambda > \lambda/2$
        and remain active under Strang;
  \item off-support indices with $\lambda/2 < |d_j| < \lambda$ are activated by Strang but
        not by LT, contributing additional positive terms to $C_D$;
  \item off-support indices with $|d_j| \leq \lambda/2$ remain inactive under both schemes.
\end{itemize}
The threshold $\lambda/2$ corresponds to the individual step size in the double proximal step
of Strang-ISTA, each of which applies $\operatorname{prox}_{(h\lambda/2)\|\cdot\|_1}$.
\end{conjecture}

\begin{remark}[Numerical evidence for Conjecture~\ref{conj:strang_support}]
\label{rem:conj_evidence}
For the synthetic test problem ($n{=}200$, $d{=}20$, $\lambda{=}0.05$,
$h \in [0.1/L,\, 0.7/L]$), the dual variable at the LT fixed point gives
$|d_8| \approx 0.563\lambda > \lambda/2$; component $j{=}8$ is activated by Strang-ISTA for
\emph{all} tested step sizes.  All remaining off-support components satisfy $|d_j| < \lambda/2$
and remain inactive.  The activated off-support component accounts for the approximately $20\%$
gap between the analytic lower bound $C_D^{\mathrm{lb}} = \|A\vect{s}^*\|^2/8 \approx 0.769$
and the empirical value $C_D \approx 0.959$.
\end{remark}

\begin{corollary}[Splitting order matters for sparse optimisation]
\label{cor:splitting_order}
Under the LASSO setting~\eqref{eq:lasso}, the LT (first-order) and Strang (second-order)
splitting schemes achieve the same $O(1/(hK))$ transient rate, but differ in their
asymptotic regime:
\begin{itemize}
  \item \textbf{LT / ISTA}: converges to $\vect{\theta}^*$ with zero asymptotic bias.
  \item \textbf{Strang-ISTA}: converges to the biased fixed point $\bar{\vect{\theta}}^h$
        with $F(\bar{\vect{\theta}}^h) - F^* = C_D h^2\lambda^2 + O(h^3)$.
\end{itemize}
Hence, for sparse optimisation problems with $\ell_1$ regularisation, the first-order
splitting (LT/ISTA) is \emph{preferable} despite the theoretical second-order accuracy
advantage of the Strang scheme for smooth ODEs.
\end{corollary}

\subsection{Numerical Experiments: LT vs Strang-ISTA for LASSO}
\label{sec:lasso_experiments}

We empirically validate Theorems~\ref{thm:lt_lasso} and~\ref{thm:strang_ista_convergence}
on a synthetic LASSO instance. We set $n = 200$, $d = 20$, $A/\sqrt{n} \in \mathbb{R}^{200\times20}$
(so $L = \|A^\top A\| \approx 1$), $\lambda = 0.05$, and a 6-sparse ground truth $\vect{\theta}^*$.
Step sizes $h \in \{0.1/L, 0.2/L, 0.3/L, 0.4/L, 0.5/L\}$ are tested; $K = 2000$ iterations.
All experiments use the estimation script \texttt{scripts/article3\_cd\_estimation.py}.

\begin{figure}[h]
  \centering
  \includegraphics[width=0.85\linewidth]{lasso_convergence.pdf}
  \caption{Convergence of LT (ISTA) and Strang-ISTA for LASSO ($h = 0.3/L$, $K = 2000$).
    LT achieves $F(\bar{\vect{\theta}}_K) - F^* \to 0$ (zero bias, Theorem~\ref{thm:lt_lasso});
    Strang-ISTA stagnates at a positive bias floor $C_D h^2\lambda^2$
    (Theorem~\ref{thm:strang_ista_convergence}).}
  \label{fig:lasso_convergence}
\end{figure}

\begin{figure}[h]
  \centering
  \includegraphics[width=0.85\linewidth]{lasso_bias_floor.pdf}
  \caption{Asymptotic bias floor of Strang-ISTA vs $h^2$.
    The linear fit (dashed) gives slope $C_D \cdot \lambda^2 \approx 0.96 \times (0.05)^2 = 2.40\times10^{-3}$
    (numerical estimate $C_D \approx 0.96$, Remark~\ref{rem:cd_formula}).
    The on-support lower bound $\|A\operatorname{sign}(\vect{\theta}^*)\|^2/8 \approx 0.769$
    (eq.~\eqref{eq:cd_formula}) underestimates $C_D$ because off-support components
    are activated at the Strang fixed point.
    LT bias floor is zero across all step sizes.}
  \label{fig:lasso_bias_floor}
\end{figure}

\begin{figure}[h]
  \centering
  \includegraphics[width=0.85\linewidth]{lasso_fixed_point_shift.pdf}
  \caption{Component-wise fixed-point shift $\bar{\vect{\theta}}^h - \vect{\theta}^*$ for Strang-ISTA
    vs the theoretical prediction $(h\lambda/2)\operatorname{sign}(\vect{\theta}^*)$
    (eq.~\eqref{eq:strang_shift}).
    Agreement is within numerical precision ($h = 0.3/L$).
    LT fixed point coincides with $\vect{\theta}^*$ (zero shift).}
  \label{fig:lasso_fp_shift}
\end{figure}

Figure~\ref{fig:lasso_convergence} confirms the qualitative prediction: LT converges
to the true optimum while Strang-ISTA stagnates at a non-zero floor.
Figure~\ref{fig:lasso_bias_floor} verifies the $h^2$ scaling of the bias.
The numerically estimated constant $C_D \approx 0.96$ exceeds the on-support lower bound
$\|A\operatorname{sign}(\vect{\theta}^*)\|^2/8 \approx 0.769$ (Remark~\ref{rem:cd_formula}),
consistent with the activation of off-support components at the Strang fixed point.
Figure~\ref{fig:lasso_fp_shift} directly confirms the fixed-point shift formula
$(h\lambda/2)\operatorname{sign}(\vect{\theta}^*)$ at the component level.

\begin{figure}[h]
  \centering
  \includegraphics[width=0.85\linewidth]{lasso_dual_variables.pdf}
  \caption{Dual variables $|d_j| = |[A^\top(A\vect{\theta}^* - \vect{y})]_j|/\lambda$ at the
    LT fixed point for all $d{=}20$ components, sorted by magnitude.
    The horizontal dashed lines mark $|d_j|/\lambda = 1$ (LT activation threshold) and
    $|d_j|/\lambda = 0.5$ (Strang activation threshold, Conjecture~\ref{conj:strang_support}).
    Support components $j \in \{0,\ldots,5\}$ achieve $|d_j|{=}\lambda$; off-support component
    $j{=}8$ has $|d_8|{\approx}0.563\lambda > \lambda/2$ and is the only additional component
    activated by Strang-ISTA.}
  \label{fig:lasso_dual_vars}
\end{figure}

\begin{figure}[h]
  \centering
  \includegraphics[width=0.85\linewidth]{lasso_support_expansion.pdf}
  \caption{Support of LT and Strang-ISTA fixed points as a function of step size $h/L$.
    LT (blue) maintains the true 6-sparse support for all tested $h$.
    Strang-ISTA (orange) activates one additional off-support component ($j{=}8$,
    $|d_8|{\approx}0.563\lambda$) for all $h \in [0.1/L,\, 0.7/L]$,
    consistent with Conjecture~\ref{conj:strang_support}.}
  \label{fig:lasso_support_exp}
\end{figure}

Figures~\ref{fig:lasso_dual_vars} and~\ref{fig:lasso_support_exp} provide evidence for
Conjecture~\ref{conj:strang_support}.  The dual variable plot (Figure~\ref{fig:lasso_dual_vars})
shows that exactly one off-support component falls in the interval $(\lambda/2,\lambda)$,
and the support expansion plot (Figure~\ref{fig:lasso_support_exp}) confirms this component is
activated by Strang-ISTA independently of step size.

\end{document}