\section{Introduction}
\IEEEPARstart{T}{he}  classification of static images is one of the central problem in Computer Vision in Machine Learning \cite{Duda2000,Cipolla2013} and has a broad range of applications such as in security and authentication with automatic facial recognition \cite{Liu2002,Ahonen2006} or object detection \cite{Ramanan2010}, in health sciences with the automatic analysis or segmentation of medical images \cite{Hames1992}, or intelligent character recognition for handwritten interpretation, such as automatic zip-code identification \cite{Lecun1990}.

Over the last decade, the level of accuracy of automatic image classification techniques has undergone a dramatic improvement coincidentally with to the breakthrough of deep-learning  and a particular class of neuro-inspired algorithms, known as convolutional neural networks (CNNs) \cite{lecun2015deep}. CNNs are multi-hidden-layered feed-forward neural network trained to extract hierarchical features from images and improve the recognition procedure thanks to their robustness to variability in viewpoint, scale, illumination conditions, and intra-class variations intrinsically present in images.

Deep learning approaches based on CNNs have outclassed simpler classification systems, such as support vector machines (SVM) \cite{cortes1995support}. As a result, in $2015$, a CNN comprised of over $150$ layers has surpassed for the first time the accuracy of humans in classifying images from the ImageNet dataset \cite{ILSVRC15}, which is made of approximately 14 millions images distributed over more than $21,000$ different categories : It reached a classification error rate of $3.57\%$ to be compared to $5.1\%$ for humans. On smaller image datasets, such as the MNIST database made of 70,000 images of handwritten digits distributed over 10 different classes, deep-learning approaches reaches a classification error as low as $0.21\%$ \cite{lecun2015deep}. However, achieving high accuracy requires the training of a large number of hidden layers (\textit{i.e.} the fine tuning of millions of parameters). This is only possible by leveraging (very) large pre-labeled databases and the high-speed, massively parallelizable, computing power from graphical processing units (GPU), thus making the training procedure computationally demanding, usually intractable on regular computers, and energy intensive.

 An alternative approach to deep-learning for image classification uses (i) simple feature extraction techniques in combination with (ii) a recurrent neural networks (RNN) and (iii) a simplified training procedure with a limited number of learnable parameters. Reservoir Computing (RC) is a Machine Learning paradigm, which includes echo-state networks (ESN) and liquid-state machines (LSM), to design and train artificial recurrent neural networks \cite{jaeger2004harnessing,maass2002real}. Reservoir Computing exploits the transient response of a randomly-, sparsely-interconnected RNN made of nonlinear dynamical discrete-time nodes, which are subjected to randomly-weighted input signals. The output layer is made of (possibly) several outputs that realize linear combinations of the nodes' transient states. The training procedure for the Reservoir Computer is limited to the optimization of the output layer weights, which is usually obtained by simple linear (or ridge) regression \cite{lukosevicius2009reservoir}. This leads to a significant alleviation in the computational complexity of the training, because only the final unidirectional connections are optimized by comparison to the entire RNN in typical training procedures for artificial neural networks. Additionally, the reduced number of trainable parameters makes the reservoir computer efficient for classification even on small databases.

 Because of its structural simplicity, reservoir computers have motivated many hardware implementations ranging from electronics \cite{appeltant2011information,Haynes2014} and spintronics \cite{Grollier2017,Grollier2019}, to photonics \cite{larger2012photonic,brunner2013parallel} with the goal of reaching ultra-fast processing speed with an energy consumption at least two order of magnitudes below that of a software-based RC running on a computer. Unprecedented classification speed have been recently achieved on the spoken-digit recognition task from the TIMIT46 data base with real-time processing speed ranging from 300,000 to 1,000,000 words analyzed per second using laser with optical feedback \cite{brunner2013parallel} and optoelectronic oscillators \cite{larger2017high}, respectively.

 The challenge in photonics-based implementations of reservoir computing stems from the experimental difficulty to couple optically a large number (typically greater than $100$) of photonic devices together. To overcome this, the time-delay reservoir computer was introduced first with electronics systems \cite{appeltant2011information} and then with various optical/optoelectronic feedback configurations \cite{brunner2013parallel,larger2012photonic,larger2017high,paquot2012optoelectronic,Uchida2018,duport2012all} with various applications notably in optical communications \cite{Argyris2018}. The principle relies on a single nonlinear dynamical node subjected to time-delay feedback, where virtual nodes are located. This approach allows theoretically to increase the size of the network linearly with the length of the feedback loop but at the expense of the processing speed decreasing because data cannot be fed to the systems faster than the time-delay to get all the virtual nodes subjected to the data \cite{Soriano2017}.

To overcome this inherent drawback of the time-delay approach, true spatiotemporal photonic reservoir computers have been proposed : (i) using photonic integrated circuits (PICs) but limited in scale to a few tens of nodes due to high losses \cite{vandoorne2014experimental,Katumba2018,katumba2019neuromorphic}, which is inadequate for tackling complex tasks such as image classification or (ii) using free-space optics with a number of nodes reaching up to 2,500 in the most recent experiments \cite{bueno2018reinforcement}, thus making possible for such architecture to solve more complex tasks.

The paper is organized as follows.
First, in Section \ref{sec:rc} we present the basics of reservoir computing, provide a description of our experimental architecture for a large-scale photonic reservoir computer implementation exceeding 10,000 neurons, and derive the corresponding physical model used in our numerical simulations. Then, in Section \ref{sec:mnist}, we introduce the particular task of handwritten digit recognition based on the MNIST database and how our reservoir computer is adapted to solve it. We also describe the various feature extraction strategies (\textit{i.e.} histograms of oriented gradients, zoning, and Gabor filters) to generate the most relevant information from the images prior to their injection in the photonic recurrent neural network.
Finally, we devote the Section \ref{sec:res} to the discussion of our results, where we compare three conceptually different modes of operation for our reservoir computer to successfully classify handwritten digits : the feedforward mode of operation, the recurrent mode with echo, and the recurrent mode with data-splitting. We analyze and compare the performance of these different modes of the reservoir computer when paired with the feature extraction layer. Finally, Section \ref{sec:ccl} concludes the paper and summarizes our contributions.

\section{Spatiotemporal Photonic Reservoir Computer}
\label{sec:rc}

In this section, we describe the basics principles of reservoir computing and how we realize and model our experimental spatiotemporal photonic reservoir computer with up to 16,384 nodes.

\subsection{Principles of reservoir computing}
\label{subsec:rc}

A reservoir computer consists of a RNN with $n\in\mathbb{N}$ discrete-time nonlinear randomly interconnected dynamical systems, which is mathematically described by the following state equation
\begin{equation}
  \mathbf{x}(k+1) = \mathbf{f}_{\text{NL}} \left( W_{res} \mathbf{x}(k) + W_{in} \mathbf{u}(k) \right),
  \label{eq:rcevo}
\end{equation}
where $\mathbf{x}(k)\in \mathbb{R}^n$ is the state vector of the RNN at discrete time $k\in\mathbb{N}$; $\mathbf{f}_{\text{NL}}$ is a nonlinear vector flow mapping the state space $\mathbb{R}^n$ to itself. $W_{res}\in\mathbb{R}^{n\times n}$ is the adjacency matrix of the RNN describing the topology of the weighted interconnections between various neurons. The spectral radius of the matrix $W_{res}$ is usually taken in the range $0.8-1.1$ to guarantee the overall stability of the RNN and the fading memory property, one of the necessary condition on the RNN to be used in the framework of reservoir computing \cite{jaeger2001echo,jaeger2004harnessing}. $W_{in}\in\mathbb{R}^{p\times n}$ is the input matrix representing the weighted interconnections between the input layer and the RNN, also known as the \textit{input mask}. Similarly to $W_{res}$, its  weights are randomly distributed with zero mean. Lastly, $\mathbf{u}(n)\in\mathbb{R}^{p}$ is the input vector.

When subjected to time-varying input data, the state vector of the RNN undergoes complex transient dynamics before returning to a stable quiescent state. These transient dynamics are wire-tapped and weighted to form an output vector $\mathbf{y}(n)$ described by the following linear relation
\begin{equation}
  \mathbf{y}(k) = W_{out} \mathbf{x}(k),
  \label{eq:rcout}
\end{equation}
where $W_{out}\in{n\times m}$ is the output matrix  representing the unidirectional interconnections between the RNN and the ouput layer. The elements of the output matrix, also called \emph{readout weights}, are the only learnable parameters of our reservoir computer architecture. They are trained using a convex optimization problem known as ridge regression (or Tikhonov regularization) \cite{tikhonov1995numerical} that reads :
\begin{equation}
   W_{out,opt} = \argmin_{W_{out}\in\mathbb{R}^{m\times n}}\|\mathbf{Y}_t-W_{out}\mathbf{X}\|^2 + \lambda\|W_{out}\|^2,

\end{equation}
where $\mathbf{Y}_{t}\in \mathbb{R}^{m\times k_{tr}}$ is the target matrix containing the $k_{tr}$ targets values for the $m$ outputs of the reservoir computer and $\mathbf{X}\in \mathbb{R}^{n\times k_{tr}}$ is the reservoir states matrix containing the RNN's transient states $\mathbf{x}(k)$ sampled during the training, when inputs vectors are fed to the reservoir during $k_{tr}$ time steps. $\lambda$ is a positive factor called ridge parameter used to control the norm of the readout weights by penalization; it is usually determined by cross-validation, when $\lambda=0$ then the ridge regression corresponds to linear regression.

The resolution of the optimization problem to obtain the optimized readout weights in $W_{out,opt}$ is performed either off-line or online \cite{antonik2017online}, \emph{i.e.} either after all the training input data are injected in the reservoir, or in the course of the injection.

\subsection{Experimental photonic reservoir computer}
\label{sec:RCexp}

The proposed photonic network is depicted in Fig{.} \ref{fig:experiment} and its structure is inspired in parts from recent experiments in Refs{.} \cite{hagerstrom2012experimental,bueno2018reinforcement}. The setup consists of a photonic arm comprising a collimated incoherent monochromatic light source at $532$ nm (Thorlabs  M530L3), a set of two polarizers oriented at a $45^\circ$ angle with respect to the vertical axis, placed before and after a liquid-crystal on Silicon (LCoS), phase-only, spatial-light modulator (SLM) with a $512\times 512$ pixels with a $8$ bit resolution (Meadowlark P$512-0532$). This particular combination of polarizing optics and a SLM allows for the phase pattern imprinted in the transverse plan of the optical beam to be nonlinearly converted into an amplitude pattern and hence ensures an all-optical implementation of the nonlinear flow $\mathbf{f}_{NL}$ from Eq{.} (\ref{eq:rcevo}). Finally, the intensity pattern is imaged through a lens on a 1.3-Million-pixels high-speed camera with 10-bit resolution (Allied Vision Mako U130B).

The detected patterns are digitized and processed in the electronic arm of the setup. First, they are linearly combined on a computer (which could potentially be replaced by a DSP board), that effectively realizes the adjacency matrix of the network $W_{res}$ with full control over their choices. Then, input signals are added numerically after the masking operation with the $W_{in}$ matrix. Finally, these signals are fed back to the SLM controller, that changes the liquid-crystal orientation of each SLM's pixel at a given discrete time step and hence generates a new phase pattern based on the signals detected at the previous time step. This setup forms a dynamical photonic-based network, which we will be described mathematically as coupled maps in Sec. \ref{sec:RCmod}.

The size of the reservoir is essentially limited by the lowest resolution between that of the SLM and the camera. In our current setup, this means that, potentially, we could reach a network size of $512\times512 = 262,144$ nodes, where each pixel of the SLM, imaged with a $1:1$ ratio on the camera, could be considered as a physical node in our photonic network. However, we only use pixels located at the center of the SLM matrix due to constrains stemming from (i) inhomogeneous intensity distribution of the optical beam and (ii) small misalignment between the SLM and camera optical axes. To mitigate the effects of these experimental imperfections, we regroup neighbouring pixels together to form bigger square-shaped pixels (see Fig. 1(b)), latter referred to as \emph{macro-pixels} in this study, and becoming our physical nodes in the photonic network. The smallest macro-pixel effectively usable in our experience is $3\times 3$, which allows for $n=16,384$ macro-pixels available at the center of the SLM. This is currently our experimental upper-bound for the network size, but using a higher-resolution camera and SLM, and possibly smaller macro-pixel sizes with more stringent alignment conditions, our network could easily scale to hundreds thousands nodes.

\begin{figure}[t!]
  \includegraphics[width=8.5cm]{figs/Figure1_JSTQE_vf.jpg}
  \caption{Scheme of the experimental setup composed of an optical arm linked to a computer. A collimated polarized green (532 nm) LED beam is reflected by the surface of the spatial light modulator (SLM). The SLM is imaged onto a camera through a polarizer (Pol{.} 2) and an imaging lens. Both the SLM and the camera are driven by a computer running a Matlab script. The latter generates the inputs from the MNIST database images, then computes the pixel values to be loaded on the SLM. Larger macro-pixels made of groups of individual SLM's pixels are used to train the reservoir in such a way to facilitate the separation on the raw camera images. The computer receives the data containing the reservoir states from the camera, computes the outputs and generates the output digits.}
  \label{fig:experiment}
\end{figure}

The speed of the setup is defined by the time needed to compute the next SLM matrix from the raw camera image. This operation is deliberately performed in Matlab for greater flexibility of the experimental scheme. The system is capable of processing 2 images per second with a larger reservoir ($n=16,384$) and up to 7 frames per second with a small reservoir ($n=1,024$).
The system's speed limitation can be alleviated by replacing the computer with a dedicated digital signal processing (DSP) board, or a field-programmable gate array (FPGA) chip, capable of performing the matrix-products computations in real time (as in e.g. \cite{antonik2017online}).
Matrix multiplication can also be offloaded to fully parallel optics \cite{bueno2018reinforcement, psaltis1985optical}.
Furthermore, a dedicated computing unit (FPGA or DSP board) could also take care of the feature extraction stage, thus allowing real-time processing of the input images.

\subsection{Physical Modelling of the photonic reservoir computer}
\label{sec:RCmod}

In this section, we propose a phenomenological modelling of the dynamics of our photonic network.

We first model the optical arm and assume, for simplicity, that the optical field $\mathbf{E}_i(t_k)$ and phase pattern $\phi_i(t_k)$ are homogeneous and constant over the surface of the $i$-th macro-pixel during $\Delta t = t_{k+1}-t_{k}$, the sampling time of our electronic feedback loop. Furthermore, using only the macro-pixels close to the optical axis of the SLM, we also assume homogeneity of the optical field over the entire $n=16,384$ macro-pixels and hence $\mathbf{E}_i(t_k) = \mathbf{E}(t_k)$ and its constant value $\mathbf{E}(t_k)=\mathbf{E}_0$ during our experiment.

The incoherent, unpolarized light source is filtered by a linear polarizer oriented at $45^\circ$ with respect to the vertical axis, hence leading to the following expression $\mathbf{E}_{0} = E_{0}/\sqrt{2}[1,1]^T$. Each SLM's  macro-pixel is composed of liquid-crystal cells on top of a highly-reflective surface. The liquid crystal cells of the $i$-th macropixel can be considered as a birefringent material with its fast-optical axis along the vertical direction and with a programmable phase-shift $\phi_i(t_k)$. Hence, the Jones Matrix for the SLM's $i$-th macropixel is given by
\begin{equation}
J_{SLM,i}(t_k) = \left(\begin{array}{cc}
     e^{i\phi_i(t_k)/2} & 0  \\
     0 & -e^{-i\phi_i(t_k)/2}
\end{array}\right).
\end{equation}
After reflexion on the SLM $i$-th macro-pixel and transmission through the second polarizer oriented at $45^\circ$ and described by the Jones matrix $J_{p,45}$, the output field is given by $\mathbf{E}_{out,i}(t_k) = J_{p,45}J_{SLM,i}(t_k)\mathbf{E}_{0}=\mathbbm{i}\frac{E_0}{\sqrt{2}}\sin\left(\frac{\phi_i(t_k)}{2}\right)\left[1,1\right]^T$.
The intensity detected by the camera and associated to the $i$-th macro-pixel hence reads
\begin{equation}
    I_i(t_k) = \|\mathbf{E}_{out,i}(t_k)\|^2 = I_0\sin^2\left(\frac{ \phi_i(t_k)}{2}\right),
\end{equation}
with $I_0 = \|\mathbf{E}_0\|^2$ the constant and uniform optical intensity emitted by the LED.

The camera has a 10 bits resolution and detects a quantified version of the intensity of the $i$-th macro-pixel on the SLM that we will denote $\lfloor{I_i(t_n)}\rfloor_{10}$. We define the $n$-dimensional state vector of the photonic network from the macro-pixel intensities detected by the camera by:
\begin{equation}
\mathbf{x}(t_k) = \left[\lfloor {I_1(t_k)} \rfloor_{10},\dots,\lfloor {I_n(t_k)}\rfloor_{10}\right]^T.
\end{equation}

The 10-bits quantified intensities detected by the camera (corresponding to the SLM's macro-pixels) are linearly combined together and masked input data is added digitally. This operation on the states are then quantified over $8$ bits before being fed back to the SLM controller so that phase values of the SLM's macro-pixels $\phi_i(t_{k+1})$ for $i={1,\dots,n}$ are updated and so are the intensities $I_i(t_{k+1})$. The output of the reservoir are determined by combining linearly the intensities detected by the cameras. This leads to the following state and output equations :
\begin{eqnarray}
 \mathbf{x}(t_{k+1}) &=& \left\lfloor I_0 \sin^2\left(\left\lfloor W_{res}\mathbf{x}(t_k) + W_{in} \mathbf{u}(t_k) \right\rfloor_{8}\right)\right\rfloor_{10},\label{eq:evoPhotRC}\\
 \mathbf{y}(t_k)&=&W_{out}\mathbf{x}(t_k).\label{eq:outPhotRC}
\end{eqnarray}
Equation (\ref{eq:evoPhotRC}) should be understood as the nonlinear function $\lfloor I_0\sin^2\left(\lfloor \cdot \rfloor_{8}\right)\rfloor_{10}$ applied to each coordinate of the $n$-dimensional phase space. Equations (\ref{eq:evoPhotRC})-(\ref{eq:outPhotRC}) are the tailored version of Eqs. (\ref{eq:rcevo})-(\ref{eq:rcout}) specific to our photonic implementation. This model will be later used in our numerical simulation to assess the performance of the photonic reservoir computer.

\section{Reservoir Computing with Feature extraction for image classification}
\label{sec:mnist}

Feature extraction is a common approach in computer vision to enhance the classification system, here the photonic reservoir computer, by providing the most relevant information in images.
In this section, we first introduce the approach adopted with reservoir computer to solve the MNIST handwritten digit recognition task and then present the five
feature extraction techniques
that we tested in this study.

\subsection{Resolution of the handwritten digit recognition task from the MNIST database with reservoir computing}

The principle of the handwritten digit recognition in the context of reservoir computing with feature extraction is illustrated in Fig{.} \ref{fig:principle}.
In this work, we used the popular MNIST database\cite{lecun1998gradient}, publicly available online, which contains $70,000$ images of handwritten digits from 0 to 9. All images have been normalized to fit into a $28\times28$ pixels bounding box, anti-aliased, and converted into gray-scale levels.

\begin{figure}[t]
  \centering
  \includegraphics[width=8.5cm]{figs/Figure2_JSTQE_vf.jpg}
  \caption{Handwritten digit recognition in the context of photonic reservoir computing with feature extraction. The images from the MNIST database undergo a feature extraction stage, where different algorithms under investigations are applied (or none, in the case of raw images). The resulting features are fed into the photonic reservoir computer with 10 binary output nodes, one for each digit. The nodes are trained to output 1 for the digit associated with the node and 0 for the other digits. The final classification is obtained by selecting the node with the maximum output, \textit{i.e.} the winner-takes-all decision strategy.}
  \label{fig:principle}
\end{figure}

The input images can be pre-processed in several ways to the input layer of the reservoir computer, depending on the type of feature extraction presented in Sub-section \ref{subsec:feat}, and it will be shown in Section~\ref{sec:res} that this choice impacts the classification error significantly.

The output layer is made of $10$ binary outputs, introduced to recognize the 10 different class of digits. Hence, each binary output is trained to give a ``1'' for an image of the digit it is associated to and ``0'' for all other digits. The winner-takes-all decision strategy is used to classify each image based on the binary output with the highest value.

Furthermore, the photonic reservoir computer can have different modes of operation and processes  the masked input data differently depending on how the latter is fed into the reservoir.
In this study, we consider three different modes of operation.

\emph{Full-image feedforward mode} -- In this mode, one full MNIST image is sent to the reservoir at each timestep $t_k$. To avoid cross-talk between consecutive images caused by the internal memory of the system, the reservoir should be exploited as memoryless map with $W_{res}=0_{n\times n}$. Therefore, the state equation of the reservoir reduces to:
\begin{equation}
    \mathbf{x}(t_{k+1})=\left\lfloor I_0 \sin^2\left(\left\lfloor W_{in} \mathbf{u}(t_k) \right\rfloor_{8}\right)\right\rfloor_{10}.
\end{equation}
The results obtained with this approach will be presented in Sec. \ref{subsec:res:ff}.

\emph{Full-image recurrent mode} -- A more advanced mode of operation is to present a full MNIST image at a given timestep $t_k$ and let the reservoir process the information during $k_e$ additional timesteps until the system relaxes into a quiescent state. Then, a new image can be sent at $t_{k+ke+1}$. In this mode of operation, $ W_{res}\ne 0_{n\times 0}$ and the state equation of the reservoir is given by Eq{.}~(\ref{eq:evoPhotRC}). We will discuss this approach further in Sec. \ref{subsec:res:rc}.

\emph{Column-wise recurrent mode} $-$ The third mode of operation for the reservoir is inspired by Ref{.}\cite{schaetti2017echo} and consists of dividing each MNIST image into $28$ columns. The image is then presented at a rate of one column per time step $k$, so that it takes 28 timesteps to feed a single image to the reservoir. In this mode of operation, the dynamics of the network is described by Eq{.}~(\ref{eq:evoPhotRC}) with $W_{res}\ne 0_{n\times 0}$. The underlying idea here is to transform the spatial correlations naturally present in images into temporal correlations in the dynamical states of the reservoir and thus exploit the internal memory of the system created by the recurrence of the network. The disadvantage, however, is a processing speed divided by $28$ compared to the first mode of operation (full-image feedforward mode). The results for this mode of operation will be presented in Section \ref{subsec:res:col}.

\subsection{Feature extraction approaches}
\label{subsec:feat}

We implemented and compared five feature extraction techniques, also considered in \cite{bahi2015robust}.
The three approaches that gave the best results are presented in this section, and their results will be discussed in Sec. \ref{sec:res}.
Table \ref{tab:feat} summarizes these techniques and their respective input dimensionalities (i.e. the number of features).
The two other approaches that did not yield acceptable classification errors are only briefly discussed in the end of this section.

\begin{table}[!htbp]
  \centering
  \begin{tabular}{lc}
     \hline
    Feature extraction method & Input dimensionality \Tstrut\Bstrut\\
    \hline
    Raw images & 784 \TstrutL\\
    Zoning 2 & 196 \\
    Zoning 4 & 49 \\
    Gabor filters & 40 -- 1152 \\
    HOG & 324 -- 1296 \BstrutL\\
        \hline
  \end{tabular}
  \vspace{0.25cm}
  \caption{Feature extraction techniques investigated in this work}
  \label{tab:feat}
\end{table}

\subsubsection{Raw images}
\label{subsec:feat:raw}

We used raw images -- \textit{i.e.} without any pre-processing -- as a benchmark for the feature extraction methods considered below.
In this particular case, grey-scale values of the pixels are used as inputs to the reservoir computer, and the dimensionality of the input is $28 \times 28 = 784$. An example of raw image is shown in Fig{.}~ Fig. \ref{subfig:feat:raw}.

\subsubsection{Zoning}
\label{subsec:feat:zon}

Introduced in Ref{.}~\cite{hussain1972results}, it is the simplest feature extraction technique considered here. It is a statistical region-based approach, where the image is divided into smaller zones and pixel densities are computed in each one. In other terms, this method consists of a combination between a convolution of the image by a filter, followed by pooling a single value from each zone of the image.

Here, we consider two variants of this approach, designed specifically for the MNIST database: \emph{Zoning 2} with $2\times2$-pixel zones, and \emph{Zoning 4} with larger $4\times4$-pixel zones.
The filters are $2\times2$ or $4\times4$ matrices filled with $1$s, thus giving a non-normalized average pixel values for each zone. This technique allows to reduce the input-dimensionality by ``dropping'' $3$ out of $4$, or $15$ out of $16$ pixels for Zoning 2 and Zoning 4, respectively. Consequently, the input dimensionality is $196$ for Zoning 2, and $49$ for Zoning 4. The resulting inputs to the reservoir computer are illustrated in Figs{.} \ref{subfig:feat:zon2} and \ref{subfig:feat:zon4}.

\subsubsection{Gabor filters}
\label{subsubsec:feat:gab}
They are linear filters mainly used for texture analysis \cite{jain1991unsupervised} and have been also used to model the receptive field profiles of simple cells in mammalian visual cortex  \cite{bovik1990multichannel,turner1986texture,perry1989segmentation,tan1990texture}. Therefore, pre-processing images with Gabor filters could be compared, to a certain extent, to perception in the human visual system.
The 2D-impulse response $h(x,y)$ of a Gabor filter is defined by the product of a Gaussian function with a sine/cosine  function. There are either complex-valued Gabor filters \cite{bovik1990multichannel} or pairs of Gabor filters with quadrature-phase relationship \cite{turner1986texture,perry1989segmentation,tan1990texture}. In this work, we use real-valued, even-symmetric Gabor filters\cite{jain1991unsupervised,malik1990preattentive}, which  impulse response is given by :
\begin{equation}
  h(x,y) = \exp \left( - \frac{x_\theta^2}{2\sigma^2_x} - \frac{y_\theta^2}{2\sigma_y^2} \right) \cos \left( \frac{2\pi}{\lambda} x_\theta \right),
  \label{eq:gabor}
\end{equation}
with
\begin{align}
    x_\theta & = x \cos (\theta) + y \sin (\theta), \\
    y_\theta & = y \cos(\theta) - x\sin(\theta),
\end{align}
where $\lambda$ is the wavelength of the cosine function with orientation $\theta$ with respect to the $x$-axis and $\sigma_{x,y}$ are the space constants of the Gaussian envelope along the $x,y$ axes, respectively.
A Gabor filter seeks for a specific frequency content in a direction $\theta$ by yielding higher values in locations visually similar to the filter itself after convolution with the image. The dimensionality of the resulting features depends on (i) the number of filters, differentiated by their lengths and orientations, and (ii) the computation of the local energy, that will be discussed in Section~\ref{subsubsec:res:ff:gab}.

\subsubsection{Histograms of oriented gradients (HOG)}
\label{subsubsec:feat:hog}

This algorithm is widely used in computer vision for object detection and locati    on \cite{bahi2015robust} in static images. It was originally proposed in Ref{.}~ \cite{dalal2005histograms} and is based on scale-invariant features transform (SIFT) descriptors \cite{lowe2004distinctive}. The HOG descriptors are obtained by creating histograms on the magnitude of images gradient $m(i,j)=(D_x(i,j)^2+D_y(i,j)^2)^{1/2}$ depending on their orientation $\theta(i,j)=\arctan(D_y(i,j)/D_x(i,j))$ with $D_x(i,j)$ and $D_y(i,j)$ the horizontal and vertical approximate gradients evaluated at pixel $(i,j)$. These gradients are computed using the two spatial filters $G_x = \left[ -1, 0, 1 \right]$ \text{and} $G_y = \left[ -1, 0, 1 \right]^T$.

A typical histogram consists of $9$ bins in the $[0^\circ,180^\circ]$ interval. It is calculated on $8\times 8$-pixel cells within the image. This ensures a dimensional reduction from $64$ gray-scale intensity values to $9$ cumulative frequencies values corresponding to the so-called HOG descriptors. To make them robust to lighting inhomogeneities, the histograms are block-normalized over larger regions of pixels within the image.
Similarly to Zoning, the number of features depends on the size of the image and the cell, the number of bins, and the block-normalisation. In the case of the MNIST database, the total number typically ranges between 234 and 1296. This will be further discussed in Sec. \ref{subsubsec:res:ff:hog}.
Figure \ref{subfig:feat:hog} illustrates the resulting gradients superimposed on top of an example image from the MNIST dataset.

\begin{figure}[!t]
  \centering
  \subfigure[]{\includegraphics[width=0.15\textwidth]{figs/mnist_raw.png}\label{subfig:feat:raw}}
  \subfigure[]{\includegraphics[width=0.15\textwidth]{figs/mnist_zon_196.png}\label{subfig:feat:zon2}}
  \subfigure[]{\includegraphics[width=0.15\textwidth]{figs/mnist_zon_49.png}\label{subfig:feat:zon4}}
  \subfigure[]{\includegraphics[width=0.15\textwidth]{figs/gabor_ex.png}\label{subfig:feat:gabor}}
  \subfigure[]{\includegraphics[width=0.15\textwidth]{figs/mnist_hog.png}\label{subfig:feat:hog}}
  \caption{Illustration of \textbf{(a)} raw MNIST image, and the feature extraction methods, applied to it: \textbf{(b)} Zoning 2, \textbf{(c)} Zoning 4, \textbf{(d)} Gabor with $\lambda=3$ and $\theta=22.5^\circ$ orientation capturing the middle stroke of the digit, and \textbf{(e)} HOG features (red) computed with $7\times7$ cells and 9 orientation bins, superimposed on top of the source image.}
  \label{fig:feat:zon}
\end{figure}

\subsection{Other methods}
\label{subsec:feat:misc}

The following two feature extraction methods were investigated in preliminary simulations, but were not included in the results (Sec. \ref{sec:res}) since they performed worse than the raw images. Moreover, both methods require the distinction between foreground and background pixels, and are thus more complex to implement on real-life images.

\subsubsection{Projection histograms}

Projection histograms (vertical and horizontal) were considered in \cite{bahi2015robust}. This approach consists in counting the number of white foreground pixels in each column or row, respectively. We used both histograms simultaneously -- that is, 56 features for each image -- and obtained classification errors of order of $12\%$.

\subsubsection{Distance profiles}

This descriptor calculates the pixel-distance between the border of the image and the first pixel of the foreground (i.e. the digit)\cite{siddharth2011handwritten}. Based on the four borders of the image, four profiles can be considered: left, right, top, and bottom. In this study, we investigated the combination of all four profiles, since it yields a descriptor with the most information on the image, encoded into 112 features. This approach works better than the projection histograms above, but we could not reduce the classification error below circa $10\%$.

\section{Results}
\label{sec:res}

The results of this study are divided into three sections, based on how the MNIST images are processed by the reservoir computer (see Sec. \ref{sec:mnist}).

\subsection{Full-image feedforward mode}
\label{subsec:res:ff}

Presenting one full MNIST image per time step $k$ is the simplest approach considered here. Besides the simplicity, it has the advantage of (1) being compatible with feature extraction techniques, presented in Sec. \ref{subsec:feat}, and (2) not increasing the number of processing time steps, since it is equal to the size of the MNIST database. The downside, however, is that it does not exploit the temporal dynamics of the reservoir computer. In fact, since the individual MNIST images are independent from each other, the inputs to the reservoir lack the temporal dependence required to make use of the recurrence of the network. Intuitively, one expects the system to perform better without memory, so that the classification of the current image would not be influenced by the previous input. This assumption was confirmed in both experiments and numerical simulations. Therefore, the recurrent reservoir computer is reduced to a simple feedforward network in this case.

Figure \ref{fig:ff} presents a summary of the numerical results obtained with the different feature extraction techniques introduced in Sec. \ref{subsec:feat}.
We simulated reservoir sizes from 1,024 up to 16,384 nodes. In all cases, the classification error decreases with the number of neurons.
At first glance, the HOG features produces the lowest classification errors, while raw images and Zoning 4 perform the worst.
We will further discuss each feature extraction methods in the following sections and present the experimental results for comparison.

\begin{figure}[b!]
  \centering
  \includegraphics[width=0.5\textwidth]{figs/fig_ff.pdf}
  \caption{Numerical results obtained with different feature extraction methods and various reservoir sizes. Raw images (unprocessed pixels) are shown for comparison. Zoning 2 technique performs slightly better, while Zoning 4 yields higher classification errors with larger reservoirs. Gabor filters significantly improve performance in comparison to other techniques, but fall short of the HOG algorithm that produced the best results in this study: $0.80\%$ with the largest reservoir of $n=16,384$ nodes.}
  \label{fig:ff}
\end{figure}

\subsubsection{Raw images}
\label{subsubsec:res:ff:raw}

Figure \ref{fig:raw} presents the numerical (green) and experimental (black) results obtained with raw images.
That is, unprocessed pixels were used as inputs to the reservoir computer.
We used these results as benchmark for other feature extraction techniques.
In simulations, raw images were classified with a $6.4\%$ error with the smallest reservoir ($n=1,024$ nodes) and $2.1\%$ with the largest reservoir ($n=16,384$ nodes).
In experiments, we obtained a $8.31\%$ and $2.85\%$ error with the smallest and the largest reservoir, respectively.

\begin{figure}[t]
  \centering
  \includegraphics[width=0.5\textwidth]{figs/fig_raw.pdf}
  \caption{Numerical (green curve) and experimental (black markers) results obtained with the raw images, without any pre-processing or feature extraction. While the experimental results do not match the simulations exactly, they present a similar trend.}
  \label{fig:raw}
\end{figure}

\subsubsection{Zoning}
\label{subsubsec:res:ff:zon}

Figure \ref{fig:ff} presents the numerical results obtained with both zoning methods -- Zoning 2 with $2\times 2$ windows and Zoning 4 with $4\times 4$ windows.

This technique, as presented in section \ref{subsec:feat}, is the simplest approach considered here, that consists in averaging pixel values over windows of fixed size. The Zoning 2 technique allows to slightly improve the results in comparison with the raw images. Since its effect can be seen as dimensionality reduction, this result indicates that the default $28\times 28$ MNIST images contain some undesirable information that can be removed to improve performance. Zoning 4 technique, on the other hand, performs worse than the raw images, which indicates that too much information is lost when averaging over $4\times 4$ windows.

Comparing the different curves in Fig. \ref{fig:ff}, one notes that both Zoning techniques provide the lowest classification error improvement in comparison with the other feature extraction techniques. Therefore, we conclude that those two techniques are not suitable for handwritten digit recognition task with reservoir computing. However, they remain of interest as they bring out characteristic spatial scales of the MNIST images.

\subsubsection{Gabor filters}
\label{subsubsec:res:ff:gab}

The case of the Gabor filters perfectly illustrates the main difficulty of the feature extraction approach to classification -- the user has to find the best set of filters for the task at hand, which often requires good understanding of the task under investigation.

At first, we tried the simpliest approach that consists in computing the energy of the convolution of the filter with the image. In this way, one obtain a single number for each filter. Preliminary simulations with 40 Gabor filters of 8 different orientations and 5 differents lengths have shown a classification error of circa 19\%, which is significantly worse than with raw images (see Fig. \ref{fig:ff}). This result suggests that global features extracted from MNIST images do not contain enough information to distinguish the digits. An intuitive example illustrating the inefficiency of global features is the particular case of separating the digits ``6'' and ``9''. Since the former is very similar to the latter after a rotation by $\pi$, the energy computed from the Gabor filters is approximately the same for both digits, making them indistinguishable from the global features point of view. Others digits, such as ``8'' and ``0'', with comparable curves, also yield energies confusing for the classifier. In order to distinguish these digits, the classifiers needs to know where exactly a particular direction is present in the image, and this is where one needs to use the local features.

At the next stage, we computed local energies by dividing images into smaller windows. This adds another variable parameter to the feature extraction technique : the size of the windows. A sensible choice requires the knowledge of characteristic scales of the MNIST images, and a reasonable first guess can be based on the insights obtained with the zoning techniques, presented in Sec. \ref{subsubsec:res:ff:zon}. Therefore, we computed the energy locally on a $7\times 7$ grid, that is, within $4\times 4$ windows. We fixed 8 directions $\theta$ -- from $0^\circ$ to $157.5^\circ$ by $22.5^\circ$ steps -- and tried different lengths of the filter -- from $2$ to $8$. The lowest classification error of $3.19\%$ was obtained with the filter length $\lambda=7$ pixels, and the smallest reservoir ($n=1024$ nodes). Local energy computation increases the number of features by a factor of $49$ in this case: that is, with 8 Gabor filters, one obtains $392$ features.

At this point, the performance of the reservoir computer is significantly better than with raw images, but the HOG approach remains unmatched (see Fig. \ref{fig:ff}). Therefore, we tried to improve the features obtained with the Gabor filters by taking inspiration from the HOG algorithm. Dalal and Triggs emphasize the importance of block-normalization: without it performance drop by $27\%$ in comparison \cite{dalal2005histograms}. In the next series of numerical simulations, we added $2\times 2$ L2-norm block-normalization. The immediate effect of this operation is the increase of the number of features by a factor of 4 -- from 392 to 1152. Remarkably, we obtain a $26.5\%$-improvement with $8$ Gabor filters of length $\lambda=3$. Simulating the system with filters of different lengths, we obtain the lowest classification error of $2.67\%$ with filters of length $\lambda=5$, and the smallest reservoir ($n=1024$ nodes). This corresponds to a $22.1\%$ improvement gained by adding the block-normalization.

In summary, Gabor filters could only achieve a circa $19\%$ classification error in the preliminary tests. However,  hand-tailoring the filters to the MNIST database (\textit{i.e.} adding local-energy computation, block-normalization, and setting the filter length to $5$), we managed to reduce the error down to $2.67\%$ with the smallest reservoir ($n=1024$) and $1.21\%$ with the largest reservoir ($n=16384$). We did not manage to outperform the HOG algorithm with Gabor filters in this study, but gained additional insights on the MNIST database and how the components of the HOG algorithm contribute to its performance (in particular, the block-normalisation).

\subsubsection{Histograms of oriented gradients}
\label{subsubsec:res:ff:hog}

The histograms of oriented gradients is a powerful feature descriptor in computer designed for the purpose of object detection \cite{dalal2005histograms}. Similarly to the Gabor filters, the algorithm requires the tuning of a few parameters for better performance. As described in Sec. \ref{subsubsec:feat:hog}, the key parameters are the cell size, the number of bins, and the size of the normalisation blocks. In the case of the MNIST classification task, the generally recommended values for the last two parameters -- 9 bins and $2\times2$ block-normalisation -- produce the best results. As for the cell size, the intuitive guess that smaller cells capture a better representation of the image does not hold here. In fact, the Zoning approach (see Sec. \ref{subsubsec:res:ff:zon}) demonstrated that more information does not always result in better classification.
Preliminary numerical simulations with $4\times4$-pixel cells produced a classification error of $2.1\%$ with the smallest reservoir ($n=1024$). Upon increasing the cell size up to $7\times7$ pixels, we reduced the error down to $1.59\%$ with the same reservoir size. Furthermore, the number of features was cut from $1296$ down to $324$. Further decrease of the number of features by limiting the number of bins down to $4$ (instead of $9$) did not improve performance.

Figure \ref{fig:hog} presents the numerical (red) and experimental (black) results obtained with the HOG features and different reservoir sizes.
In simulations, the MNIST digits were classified with a $1.59\%$ error with the smallest reservoir ($n=1,024$ nodes) and $0.8\%$ with the largest reservoir ($n=16,384$ nodes).
In experiments, we obtained a $2.53\%$ and $1.03\%$ error with the smallest and the largest reservoir, respectively.

These results are of the same order of magnitude as the best performance found in the literature, summarized in Tab. \ref{tab:sota}. At the moment of writing this article, the best performance on the MNIST task was obtained with CNNs \cite{wan2013regularization} with numerous approaches reporting  classification errors below $1\%$ (they will not be exhaustively presented here). We selected and listed several alternative methods to CNN, not relying on neural networks, such as SVM, k-means, k-nearest neighbours (KNN), and boosted stumps. Table \ref{tab:sota} shows that our photonic RC outperforms, in numerical simulation, several deep approaches, such as deep Boltzman machines (DBM). Furthermore, the experiment produces lower error than the earlier versions of CNNs, such as LeNet-4\cite{lecun1998gradient}.
One should keep in mind the difference between a shallow neural network, implemented in photonic hardware, and a complex convolutional neural network, implemented on ideal noiseless digital processor. Finally, the reservoir computer significantly improves the performance compared to non-linear and linear classifiers.

\begin{figure}[t]
  \centering
  \includegraphics[width=0.5\textwidth]{figs/fig_hog.pdf}
  \caption{Numerical (red curve) and experimental (black markers) results obtained with the histograms of oriented gradients.
  Similarly to the raw images, one observes a certain offset with the small and the midsize reservoir. The numerical and experimental performances obtained with the largest reservoir are fairly close, on the other hand.}
  \label{fig:hog}
\end{figure}

\begin{table}
  \centering
  \begin{tabular}{lll}
    \hline
    Authors & Method & Performance \Tstrut\Bstrut\\
    \hline
    Wan et al. \cite{wan2013regularization} & CNN & $0.21\%$ \TstrutL\\
    Wang et al. \cite{wang2016unsupervised} & K-Means & $0.35\%$ \\
    Keysers et al. \cite{keysers2007deformation} & K-NN & $0.52\%$ \\
    Decoste et al. \cite{decoste2002training} & SVM & $0.56\%$ \\
    \textbf{This work (simulations)} & \textbf{Photonic RC} & \textbf{0.8\%} \\
    K\'egl et al. \cite{kegl2009boosting} & Boosted Stumps & $0.87\%$ \\
    Schaetti et al. \cite{schaetti2017echo} & RC & $0.93\%$ \\
    Salakhutdinov et al. \cite{salakhutdinov2009deep} & DBM & $0.95\%$ \\
    \textbf{This work (experiment)} & \textbf{Photonic RC} & \textbf{1.03\%} \\
    LeCun et al. \cite{lecun1998gradient} & LeNet-4 & $1.1\%$ \\
    LeCun et al. \cite{lecun1998gradient} & Non-linear classifier & $3.3\%$ \\
    LeCun et al. \cite{lecun1998gradient} & Linear classifier & $7.6\%$ \BstrutL\\
    \hline
  \end{tabular}
  \vspace{0.25cm}
  \caption{Performance of various state-of-the-art digital approaches compared to our best results on the MNIST database (with $n=16,384$).}
  \label{tab:sota}
\end{table}

\subsection{Full-image recurrent network}
\label{subsec:res:rc}

This second approach is similar to the first one (see Sec. \ref{subsec:res:ff}) in the way that the reservoir computer is presented with the full MNIST image and the features extracted from it. The difference lies in the fact that each subsequent image is input with a fixed delay, so that the temporal dynamics of the reservoir can be exploited to process the information. The transients, i.e. the states of the nodes induced by the input, are recorded and concatenated together into a single reservoir state, used for training. This approach effectively increases the number of trainable parameters, since the size of $W_\text{out}$ is now proportional to the duration of the transients. Consequently, it also lengthens the overall processing and makes the training process more computationally expensive, especially the matrix inversion operation.

The main downside of this technique is the increased number of hyper-parameters of the reservoir. That is, on top of the input gain, the properties of the $W_\text{res}$ matrix (such as the spectral radius, the amplitude of the diagonal elements, and others) have to be optimized to improve the results. As experimental optimization of the hyper-parameters requires long runs in practice, we only investigated this method in numerical simulations.

Figure \ref{fig:rc_all} presents the numerical results obtained with this approach. The classification error is plotted against the number of transients recorded for training using the same feature extraction techniques as in Sec. \ref{subsec:res:ff}. All results were obtained with the smallest reservoir ($n=1024$) and optimized hyper-parameters. The left-most point of each curve corresponds to the feedforward network approach, that is, no transient is recorded and the images are input at each time step $k$. The graph reflects the two major observations produced by this study: (1) the additional transients improve the classification performance for all feature extraction techniques (including raw images), and (2) two transients is the optimal setting for all cases. Table \ref{tab:rc_gains} contains the performance improvements obtained with 2 transients as compared to the feedforward network.

\begin{figure}[t]
  \centering
  \includegraphics[width=0.5\textwidth]{figs/fig_rc.pdf}
  \caption{Numerical results obtained with the recurrent network ($N=1024$) approach and various feature extraction methods. Different lengths of transients have been investigated: zero transients corresponds to the feedforward network case, discussed in Sec. \ref{subsec:res:ff}. All curves decrease until two transients and then remain constant.}
  \label{fig:rc_all}
\end{figure}

\begin{table}
  \centering
  \begin{tabular}{lccc}
    \hline
    Feature extraction method & Feedforward & Recurrent & Gain \Tstrut\Bstrut\\
    \hline
    Raw images & $0.0643$ & $0.0539$ & $16.2\%$ \TstrutL\\
    Zoning 2 & $0.0500$ & $0.0417$ & $16.6\%$ \\
    Zoning 4 & $0.0547$ & $0.0455$ & $16.8\%$ \\
    Gabor filters & $0.0264$ & $0.0198$ & $25.0\%$ \\
    HOG & $0.0159$ & $0.0133$ & $16.4\%$ \BstrutL\\
    \hline
  \end{tabular}
  \vspace{0.25cm}
  \caption{Performance improvement offered by the recurrent approach with the smallest reservoir ($N=1024$).}
  \label{tab:rc_gains}
\end{table}

Since the additional transients improve the classification error for the smallest reservoir ($N=1024$), we investigated its effects on larger reservoirs. For simplicity, we only considered the HOG features, since they yield the lowest error rate with the feedforward network (see Sec. \ref{subsubsec:res:ff:hog}). At this stage, we could perceive the increased complexity of training a reservoir with transients, and could only perform numerical simulations up to $N=6400$ nodes on our desktop computers with 32 Gb of RAM. Figure \ref{fig:rc_hog} compares the classification errors obtained with feedforward and recurrent networks. The optimal number of transients was found to be 2, similarly to the previous observations. Table \ref{tab:rc_hog_gains} lists the performance gains for different reservoir sizes.

\begin{figure}[t]
  \centering
  \includegraphics[width=0.5\textwidth]{figs/fig_rc_hog.pdf}
  \caption{Numerical results obtained with the histograms of oriented gradients, processed by a feedforward network (blue) and a recurrent network (green) with 2 transients. Due to computational complexity, we could only train reservoirs up to $N=6400$ nodes.}
  \label{fig:rc_hog}
\end{figure}

\begin{table}
  \centering
  \begin{tabular}{lccc}
  \hline
    Reservoir size ($N$) & Feedforward & Recurrent & Gain \Tstrut\Bstrut\\
    \hline
    1024 & $0.0159$ & $0.0133$ & $16.4\%$ \TstrutL\\
    2304 & $0.0127$ & $0.0113$ & $11.0\%$ \\
    4096 & $0.0112$ & $0.0096$ & $14.3\%$ \\
    6400 & $0.0089$ & $0.0085$ & \hfill $4.5\%$ \BstrutL\\
    \hline
  \end{tabular}
  \vspace{0.25cm}
  \caption{Performance improvement offered by the recurrent approach with the histograms of oriented gradients and different reservoir sizes.}
  \label{tab:rc_hog_gains}
\end{table}

In summary, although the recurrent approach increases the number of trainable parameters, lengthens the training process and makes it more computationally expensive, it offers interesting results that have to be taken into account for improving the performance of our reservoir.

\subsection{Column-wise recurrent mode}
\label{subsec:res:col}

The third approach, inspired by \cite{schaetti2017echo}, uses an explicit temporal encoding of the spatial visual information from the MNIST images in order to activate the recurrent dynamics of the reservoir computer. The idea here is to split the full image into smaller portions and feed them sequentially into the classifier. The separation can be done in various ways: columns or rows (overlapping or adjacent), sliding windows, sub-images, among others. Similarly to \cite{schaetti2017echo}, we consider non-overlapping columns, which transforms each image into 28 inputs of 28 dimensions. Therefore, such temporal encoding reduces the input dimensionality, but increases the processing time proportionally.

After defining the encoding procedure, the specifics of the reservoir training have to be considered. That is, as each image corresponds to 28 inputs, which in turn give birth to 28 reservoir states, one needs to define how to use those states to train the output weights. We chose to investigate two approaches here:
\begin{itemize}
    \item The first training procedure consists in considering each column as an individual input and train the reservoir computer to output the target digit on every input column forming the corresponding image. Intuitively, one does not expect the system to produce the correct class at the first input column, but assumes the output to converge to the target digit with the last columns of the image. Therefore, while all 28 reservoir states are used for training for each image, the evaluation should only be performed on the basis of the last columns. In this work, we chose to select the last 7 columns and assign the final decision to the most frequent class among the 7 outputs.
    \item The second training method is inspired by \cite{schaetti2017echo} and consists in combining a selection of reservoir states at specific time steps into one, and computing the output weights based on the aggregate reservoir state. This idea is fairly similar to the recurrent network approach (see Sec. \ref{subsec:res:rc}), as it effectively increases the network size by combining several reservoir states into one, and thus expanding the number of trainable parameters.
\end{itemize}

The results of these two approaches are presented in Secs. \ref{subsubsec:res:col:each} and \ref{subsubsec:res:col:comb}, respectively. We used a small reservoir computer ($n=1024$) in all simulation for practical reasons.

\subsubsection{Training with individual reservoir states}
\label{subsubsec:res:col:each}

The first training option appears more challenging for the reservoir computer, since it requires a system with memory long enough to keep all relevant columns of the image.
In fact, we obtained a relatively high classification error of $21.2\%$. Upon investigation of the reservoir signals we noticed that while the outputs converge, as expected, they often converge to a wrong class. This observation suggests that individual columns do not contain enough relevant information to allow the classifier to make the correct decision.

\subsubsection{Training with aggregation of reservoir states}
\label{subsubsec:res:col:comb}

The second training option alleviates the issue encountered with single columns, as the reservoir computer is trained on an aggregation of several reservoir states activated by their corresponding input columns. This methods also relaxes the memory requirements, as the explicit aggregation of reservoir states no longer forces the network to keep the input information in memory.

The choice of the reservoir states to aggregate for the training is crucial in this approach. That is, training the system on the combination of the first three columns, for instance, is unlikely to produce low classification errors. For this reasons, we started by fixing the hyperparameters of the reservoir computer to reasonably good values (yet not the most optimal), and tested, in numerical simulations, all possible combinations of training with two reservoir states (that is, 378 trials). We found that the combination of the 17-th and 24-th reservoir states yielded the lowest error. Furthermore, we observed that no states below 10-th were chosen in the top ten results. We thus conclude that the earlier reservoir states (from the first to circa the 10-th) do not contain enough information on the digit to allow precise classification, and the later states (from 11-th to the last) should be privileged.
Then, we performed the same series of trials with 3 reservoir states.
In order to reduce the number of trials from 3276 down to 969, we ignored the reservoir states from 1 to 10, based on the observation above.
We found the best combination of states to be the 14-th, the 18-th, and the 26-th.
Finally, we tried the combinations of 4 states, and found the optimal one to be the 14-th, the 16-th, the 20-st, and the 24-th.

After defining the optimal combinations of reservoir states to aggregate for the training, we kept them fixed, and performed the optimisation of the hyperparameters. Figure \ref{fig:cols_comb} presents the classification errors with optimal training and reservoir settings. For completeness, we also evaluated the scenario with one column, and obtained a $11.7\%$ error with the 17-th column. The best performance of $4.42\%$ was produced with the combination of 4 reservoir states. Compared to the results obtained with raw images, discussed in Secs. \ref{subsubsec:res:ff:raw} and \ref{subsec:res:rc}, this method yields lower errors and thus presents the best approach on raw images among those considered in this study. Nevertheless, the $4.42\%$ classification error is no match for the HOG technique (see Sec. \ref{subsubsec:res:ff:hog}) or the state-of-the-art.

\begin{figure}[t]
  \centering
  \includegraphics[width=0.5\textwidth]{figs/fig_cols_comb.pdf}
  \caption{Numerical results obtained with column-wise input data encoding and aggregation of several reservoir states for the training process, with $n=1024$. In particular, training with $n$ corresponds to optimising the readouts weights on the basis of the sole reservoir state, produced after the processing of the 17th column of the image. Conversely, training with $4n$ consists of aggregating four reservoir states, obtained by processing the columns 14, 16, 20, and 24 of the image, for the training process.}
  \label{fig:cols_comb}
\end{figure}

\section{Conclusion}
\label{sec:ccl}

In summary, we report two major results in this paper. First, we present a large-scale experimental photonic reservoir computer, based on off-the-shelf components, capable of implementing neural networks up to $16,384$ nodes. Second, the computational power offered by our large dynamical system allows us to tackle an important task in computer vision: the recognition of handwritten digits. To this end, we study three conceptually different approaches of applying our reservoir computer to this task. First, we exploit the reservoir in a memoryless, feedforward mode of operation paired with popular handcrafted feature-extraction techniques (Zoning, Gabor filters, and HOG) prior to feeding the input layer of the photonic reservoir. The second mode of operation consisted in using the temporal dynamics of the network and hence benefit from its internal memory enabled by recurrent connections, still paired to feature extraction. Third, we used the temporal dynamics of our reservoir computer with a time-dependent signal, created by splitting the raw MNIST images into individual columns, thus creating temporal dependence and stimulating the recurrence of the system. We report a $0.8\%$ classification error in numerical simulations and $1.03\%$ error in experiments, obtained with the feedforward mode of operation using the HOG feature extraction technique. We also demonstrated in numerical simulations the benefit of the recurrence mode of operation allowing for relative gain improvement ranging from $4.5\%$ to $16.4\%$ for the smallest reservoir in our study and the superiority of preserving spatial correlation and feature extraction compared to transforming the 2D MNIST images into multiple 1D signals. In this paper, we have shown that large-scale photonic reservoir computer can perform with level of performance comparable to the best digitally-implemented neuro-inspired currently available and hence represent a meaningful  alternative to these methods, hence opening more in-depth research endeavours on exploiting photonic reservoir computer for advanced image processing.

In perspective, our work opens several directions for future research. The handcrafted features could be replaced by transfer learning, i.e. by automatically synthesised features from the first layers of a convolutional neural network trained on the MNIST dataset \cite{olivas2009handbook}. Furthermore, since the reservoir computer yields competitive results on handwritten digit recognition, one could increase the complexity of the task and apply the system to visual object recognition using e.g. the CIFAR-10 dataset \cite{krizhevsky2009learning}.
The present work is thus the first step in what we hope should be a very fruitful line of investigations.

\section*{Acknowledgment}

This work was supported by the AFOSR (grants No. FA-9550-15-1-0279 and FA-9550-17-1-0072), the R\'egion Grand-Est. The authors gratefully acknowledge Dr. Daniel Brunner for the insightful scientific exchanges and discussions related to spatiotemporal photonic reservoir computers.

\ifCLASSOPTIONcaptionsoff
  \newpage
\fi

\bibliographystyle{IEEEtran}


\section{Introduction}
XTE J1118+480 was discovered on 2000 March 29 by {\it All Sky Monitor} (ASM) onboard RXTE satellite at R.A.= $11^h18^m10^s.79$,
Dec.= $48^\circ02'12''.42$ (Remillard et al. 2000). Follow up observations were carried out by Uemura et al., 2000;
Chaty et al., 2000; Wren \& McKay, 2000; Mauche et al., 2000; Pooley \& Waldram, 2000. Because of its unique position
in the Galactic halo, this black hole binary (BHB) suffered a very less absorption and is widely studied in multi-wavelength
during both of its active outbursting periods and in the quiescence state (Revnivstev et al. 2000; Cook et al., 2000; Garcia et al., 2000;
Haswell et al., 2000; Hynes et al., 2000; McClintock et al., 2000; Wagner et al., 2000; Taranova et al., 2000; Esin et
al., 2001; Markoff et al., 2001; Hynes et al., 2003; Chaty et al., 2003; Shahbaz et al., 2005). The distance of this system
is estimated to be around 1.8 kpc (McClintock et al. 2001a). The distance above the Galactic plane ($\sim1.7~kpc$) places
it in the `Lockman halo' (Uemura et al. 2000).
Quasi Periodic Oscillations (QPOs) were observed in X-ray, optical and EUV (Wood et al. 2000; Haswell et al. 2000).
The mass of the black hole (BH) has been predicted dynamically by Wagner et al. (2001) ($6.0-7.7~M_\odot$), Gelino et al. (2006)
($8.53\pm0.6~M_\odot$), Khargharia et al. (2013) ($6.9-8.2~M_\odot$).
This is suggested to be a high inclined ($\sim68-79^\circ$) (Khargharia et al. 2013), short orbital period ($\sim4.1~hrs$)
(Patterson et al. 2000; Gonzalez Hernandez et al. 2012) low mass binary system.
The spectral and the temporal properties of the source during its 2000 outburst have been studied by Chatterjee et al. (2019; hereafter Paper-I),
using RXTE PCA and HEXTE combined data with the physical two component advective flow (TCAF) model (Chakrabarti \& Titarchuk 1995) to understand
the accretion flow dynamics of the source. Estimation of X-ray contribution of the jet made from spectral analysis with the TCAF model
suggests that the outburst was dominated by the jets/outflows. Monotonic evolution of low frequency QPOs is also studied.
In Paper-I, the probable mass of the black hole (BH) is also estimated from the spectral analysis in the range
of $6.25-7.49~M_\odot$ or $6.99^{+0.50}_{-0.74}~M_\odot$.

Roughly after five years of the quiescence phase on 2005 Jan., XTE~J1118+480 exhibited a new outbursting activity of shorter duration
($\sim 1.5$~months) and lower intensity flux compared to its earlier (2000) outburst. Multi-wavelength studies of the source
during this outburst is also carried out by many authors. Zurita et al. (2005a) reported an optical outburst on January 9, 2005 .
The outburst was also detected at X-ray and radio wavelengths (Remillard et al. 2005; Pooley 2005; Rupen et al. 2005).
The outburst faded rapidly, and reached near quiescence by the late February (Zurita et al. 2005b). The evolution of long-term
lightcurve and outburst properties had been discussed by Zurita et al. (2006).

In the present {\it paper}, our goal is to study accretion flow properties of the source during its 2005 outburst.
We wanted to check whether nature of the source is similar to other shorter orbital period harder or type-II transient
low-mass X-ray binaries recently studied by our group (see, Debnath et al. 2017 and references therein).
Here, we also make an effort to study properties of jets, particularly in X-ray band. A detailed spectral analysis is made
using physical TCAF model (Chakrabarti \& Titarchuk 1995), which is based on transonic flows (Chakrabarti, 1996).
This TCAF configuration consists of two types of flows, namely, high viscous, high angular momentum Keplerian
flow along the equatorial plane and low viscous, low angular momentum sub-Keplerian flow enveloping the Keplerian disk. Close to
the BH, axisymmetric shock forms due to the centrifugal barrier and the supersonic sub-Keplerian flow suddenly jumps to the subsonic
branch. The post-shock region puffs-up and becomes hot due to slowing down of the matter at the shock location. This puffed-up region
between the shock and the inner sonic point just outside of the horizon acts as the Compton cloud. This is known as the CENtrifugal barrier
supported BOundary Layer or CENBOL. Thermal soft photons from the Keplerian disk becomes hard via repeated inverse-Compton scattering
in the region. Although this generalized accretion flow model was introduced more than two decades ago, its recent implementation
(after generation of model {\it fits} file using $\sim 10^6$ number of theoretical spectra produced by varying five model input parameters)
as an additive table into HEASARC's spectral analysis software package XSPEC allowed one to get a more clear picture about the
flow dynamics of several BH sources (Debnath et al. 2014, 2015a,b, 2017; Mondal et al. 2014, 2016; Jana et al. 2016, 2020a;
Chatterjee et al. 2016, 2019; Molla et al. 2017). Masses of a few black hole candidates (BHCs) have been measured in a better
accuracy from normalization independent spectral analysis with the TCAF model (Molla et al. 2016, 2017; Chatterjee et al. 2016, 2019;
Jana et al. 2016, 2020a; Debnath et al. 2017; Bhattacharjee et al. 2017; Shang et al. 2019).
Theoretically, the coupling of disk and jet connection have been studied based on the transonic flow model by several authors (Chakrabarti
1999a,b; Chattopadhyay \& Das 2007; Aktar et al. 2015). The outflow rate from the inflowing accretion rate have been calculated using
hydrodynamics of infalling and outgoing transonic flows (Chakrabarti 1999a,b; Das \& Chakrabarti 1999). Chakrabarti \& Mandal (2006)
studied the two component accretion flows in the presence of synchrotron. Estimation of the contribution of the jet X-ray fluxes (if present)
and their properties are also studied with the TCAF model (Jana et al. 2017, 2020b; Paper-I). Even frequencies of the dominating
type-C QPOs are also predicted from the TCAF model fitted shock parameters (Debnath et al. 2014; Chatterjee et al. 2016).

This {\it paper} is organized in the following way: in \S 2, we discuss the observation and data analysis procedures.
In \S 3, we present the results from our spectral analysis using two types of models. Evolution of X-ray contribution of jets/outflows
during the outburst is studied with that of radio fluxes.
Finally, in \S 4, we summarize the result and present briefly a comparison between this and the 2000 outburst.

\section{Observations and Data Analysis}
In 2005, although the outburst was reported on Jan. 9 by Zurita et al. (2005a), RXTE started monitoring it
five days later on a daily basis. To find the broad spectral information in $3-100$~keV band, we studied 21 observations
of combined data of RXTE PCA and HEXTE instruments during 2005 January 14 (MJD=53384.99) to January 25 (MJD=53395.59).
HEASARC's software package HeaSoft version HEADAS 6.16 and XSPEC version 12.8 had been used for our analysis.
We followed the methods as mentioned in Debnath et al. (2013, 2015a) for data reduction and analysis.

RXTE PCA lightcurve in the energy range of $2-15$~keV and $2-25$~keV are generated using the event mode data of maximum time
resolution of $125~\mu s$. The power density spectra (PDS) are generated using XRONOS task `powspec' on 0.01 sec binned lightcurves.
We used 1 sec time binned background subtracted lightcurves of the proportional counter unit 2 (PCU2) to calculate
average PCA count rate in 2-25 keV for each observation.

In general, PCU2 data of `standard 2' mode (FS4a* in the energy range $3-25$~keV) and HEXTE science mode data (FS52* in the
energy range of $20-100$~keV) of Cluster 0 or A are used for spectral analysis. For some observations due to low s/n, we
consider HEXTE data only in $20-40$~keV band. The background subtracted spectra was first fitted with a single power-law (PL) model,
verifying that no significant thermal disk blackbody component required. After that, we refitted all the spectra using the TCAF solution
based additive table {\it fits} file. A fixed hydrogen column density $N_H=1.3\times10^{20}~atoms~cm^{-2}$ is used for photoelectric
absorption model {\it phabs} (Garcia et al. 2000; McClintock et al. 2001b). A 0.5\% systematic error is considered during spectral
analysis with the PL as well as the TCAF model.

To fit a spectrum with the TCAF model, one requires to supply four flow parameters. These parameters are: two type of accretion
rates: (i) Keplerian disk rate ($\dot m_d$ in Eddington rate $\dot M_{Edd}$), (ii) sub-Keplerian halo rate ($\dot m_h$ in
$\dot M_{Edd}$), and two shock parameters : (iii) shock location ($X_s$ in Schwarzschild radius $r_s=2GM_{BH}/c^2$),
(iv) compression ratio ($R=\rho_+/\rho_-$, where $\rho_+$ and $\rho_-$ are the densities of post- and pre- shock regions respectively).
In addition, if the mass of the black hole ($M_{BH}$ in $M_\odot$) is known, one needs to provide the value of the mass and keep as frozen.
Otherwise, each fitted spectrum returns a value of the mass. One also gets a normalization ($N$) from each fitting. This normalization
parameter is a function of the mass of the black hole ($M_{BH}$), the distance ($D$) of the system, and the inclination angle ($i$) of the disk.
So for a specific black hole, this normalization value is expected to be a constant parameter, if measured with a given instrument
(Molla et al. 2016, 2017), excluding the facts that the system is not precessing and there is no significant jet/outflow from the disk.
Jana, Chakrabarti \& Debnath (2017; hereafter JCD17) developed a method to estimate X-ray contribution of jets/outflows from
spectral analysis with the TCAF model. They compared the variation of the model normalization with the observed radio fluxes, since
radio emission is known to be a tracer of jets/outflows. If there is a significant X-ray contribution from jets, we may require
a higher value of the normalization to fit a spectrum to take care of extra photon contributions from the base of the jets.
While comparing $N$ with the radio flux ($F_R$), if lower $N$ is found when $F_R$ is also at its minimum, then it implies
that during this low normalization days there is insignificant or no X-ray contribution from the jets in the net emission spectrum.
In other words, a large amount of outflow results in higher values of $N$. They also prescribed the procedure to calculate
the amount of the X-ray outflow flux ($F_{ouf}$). Following their method in Paper-I, Chatterjee et al. (2019) calculated jet X-ray fluxes
during the 2000 outburst of XTE J1118+480, where the minimum  value of $N$ is observed at $4.36$. Since the current outburst of XTE~J1118+480
is also highly jet dominated, using the same method, we estimated the contribution of the jet X-ray flux in total flux for all observations.

Note: Since the main goal of this paper is not to measure the mass of the BH, we kept mass of the BH frozen at $7~M_\odot$ while fitting
energy spectra with the TCAF model. This is the mean value of the TCAF model fitted mass values obtained from the spectral analysis
of the 2000 outburst of XTE~J1118+480 (see, Paper-I).

\section{Results}
We study accretion properties of XTE J1118+480 during its 2005 outburst by analyzing $21$ RXTE PCA and HEXTE observations, selected
from 2005 January 14 (MJD=53384.99) to January 25 (MJD=53395.59). Both temporal and spectral properties of the source are
studied. Spectral study is done with the phenomenological (PL) and physical (TCAF) accretion flow models.

In Fig. 1, one example of power density spectrum (PDS) is shown along with the TCAF fitted combined PCA plus HEXTE spectrum of the
data with observation ID 90011-01-01-04. In Fig. 2, variation of PCU2 count rate (in $2-25$~keV) and spectral parameters fitted
with two different types of models: PL (PL flux, PL photon index $\Gamma$) and TCAF (Keplerian disk rate $\dot m_d$, sub-Keplerian
halo rate $\dot m_h$, shock location $X_s$ and compression ratio $R$) are shown.
Since we have been able to separate total observed X-ray into its two constituents: contribution from $i)$ inflowing matter or
accretion disk, and $ii)$ outflowing matter or jets; in Fig. 3(a)-(c) we show variation of total X-ray flux ($F_{X}$), X-ray fluxes
from accretion disk ($F_{inf}$) and jets ($F_{ouf}$). The variation of model normalization ($N$) and observed radio flux ($F_{R}$)
are also plotted in Fig. 3(d)-(e). The radio flux ($F_R$ in mJy) data of $15$~GHz Ryle Telescope are shown here, are
adopted from Brocksopp et al. (2010).

\begin{figure}[h]
\vskip 0.1cm
        \centerline{
        \includegraphics[scale=0.6,width=7truecm,angle=0]{fig1a.eps}\hskip 0.3cm
        \includegraphics[scale=0.6,width=7truecm,angle=0]{fig1b.eps}
        }
\caption{(a) A sample power density spectrum of $0.01$~sec time binned PCU2 lightcurve, and (b) TCAF model fitted
combined PCA plus HEXTE spectrum in $3-100$~keV energy range for observation ID 90011-01-01-04 (MJD=53387.80) are shown.}
       \label{fig1}
\end{figure}

The detailed spectral analysis results are presented in Tables 1 \& 2. In Table 1, both PL and TCAF model fitted spectral
parameters are presented. Different types of X-ray fluxes ($F_{X}$, $F_{inf}$, $F_{ouf}$), the percentage of the contribution
of the jet X-ray in total X-ray, and TCAF model fitted normalization parameter values are presented in Table 2.

\subsection{Temporal Study}
RXTE PCA started to monitor the source five days later of the report of the outburst on 2005 Jan. 9 by Zurita et al. (2005a).
Since this is a Fast Rise Slow Decay (FRSD) type outburst (see, Debnath et al. 2010), during the first PCA observation, the source
already crossed its peak flux (see, Fig. 2(a)). It showed a monotonic  decrease in the PCU2 photon count rate in $2-25$~keV energy
band. Fourier transformed PDSs (see, Fig. 1(a)) are generated for all observations using $0.01$~sec time binned $2-15$~keV lightcurves.
We have not seen any prominent signature of the low frequency QPOs, which are generally observable during the outburst of
this type of transient BHCs. Contrary to this, a monotonic evolution (increasing frequency from $0.06$ to $0.16$~Hz) of the
QPOs were observed during 2000 outburst of the source (see, Paper-I). However, there is a report of the presence of a mHz QPO
by Shahbaz et al. (2005) in the optical waveband on June 2003, which is the quiescent phase between 2000 and 2005 outbursts.

\begin{figure}[h]
\vskip 0.1cm
\centerline{
            \includegraphics[scale=0.6,width=14.0truecm,angle=0]{fig2.eps}
           }
\caption{Variation of (a) PCU2 count rate (counts/sec) in $2-25$~keV energy range, PL model fitted (b) PL flux in unit of
$10^{-9}~erg~cm^{-2}s^{-1}$, (c) PL photon index ($\Gamma$), TCAF model fitted (d) disk rate ($\dot{m}_d$) in $\dot{M}_{Edd}$
(e) halo rate ($\dot{m}_h$) in $\dot{M}_{Edd}$, (f) shock location ($X_s$) in $r_s$, (g) compression ratio ($R$) with day
(in MJD) are plotted.}
\label{fig2}
\end{figure}

\subsection{Spectral Study}
The spectral analysis is done independently with the two types of models: $i)$ phenomenological PL model, and $ii)$ physical TCAF
model. While PL model fit gives us an overview about the spectral properties of the source during different phases of the outburst,
TCAF model fitted physical flow parameters allow us to get a the accretion flow dynamics and evolution of the flow geometry.

First, we fitted the background subtracted spectra solely with the PL model. No significant disk blackbody contribution
was noticed during the spectral fittings. Although low values of the PL photon indices ($\Gamma \sim 1.75-1.87$) are
observed during the outburst, in the late declining phase, a slow increase in $\Gamma$ values (from $1.81$ to $1.87$) are
observed (see, Fig. 2(c)). On the first observation (MJD= 53384.99) day, $\Gamma=1.75$ is observed, and it varied in a
narrow range of $\sim 1.75-1.81$ for the next seven days till MJD=53392.05 before rising slowly in the next three days
to reach its maximum value of $1.87$ on the last day (MJD=53395.59) of our observation. The observation of the low $\Gamma$
signifies a hard spectral state, which is similar to what was seen during its earlier (2000) outburst of this source
(Chaty et al. 2003; Paper-I). PL flux showed a similar nature of monotonic decrease in PCA count rate during the entire
phase of the outburst (see Fig. 2(a) and 2(b)).

To obtain a physical picture of the accretion flow dynamics of the source during its 2005 outburst, we refitted all spectra
with the TCAF model based {\it fits} file as an additive table model in XSPEC. During spectral fit, we kept mass parameter
frozen at $7~M_\odot$ as this was the estimated mass value obtained from spectral analysis with the TCAF model during 2000
outburst of the source (Paper-I). Throughout the observation period, the Keplerian disk rate ($\dot m_d$) was observed to be
very less as compared to the sub-Keplerian halo rate ($\dot m_h$). This result is consistent with the  non-requirement of the
disk blackbody (DBB) component while fitting spectra with the PL model. On the first observation (MJD=53384.99) day, the value
of $\dot m_d$ was $0.041$ $\dot M_{Edd}$. On the second observation day, it increased slightly ($\sim 0.044$ $\dot M_{Edd}$),
and subsequently became roughly constant in that value except during the last few observations when a slow monotonic increase
(up to $0.48$~$\dot M_{Edd}$) was observed. The sub-Keplerian halo rate ($\dot m_h$) showed a monotonic decrease from $0.31$ to
$0.26$ $\dot M_{Edd}$ during the entire period of the outburst. The variations of $\dot m_d$ and $\dot m_h$ are plotted in
Fig. 2(d)-(e). The variation of the TCAF model fitted shock parameters (location $X_s$ and compression ratio $R$) are shown
in Fig. 2(f)-(g). On the first observation day, a strong shock ($R=3.61$) at a far location ($X_s=406~r_s$) was observed.
As the day progresses, shock becomes weaker and moves inward. On the last observation day, a comparatively weaker shock
($R=3.40$) at $X_s=245~r_s$ is observed.

\subsection{Jet X-ray}
Radio emission observed in BHCs is considered to be a tracer of jets and outflows. The strong radio emission is reported during
most part of the outburst (see, Hynes et al., 2006; Maitra et al., 2009; Brocksopp et al. 2010). During the 2000 outburst,
similar strong radio emission is observed (see, http://www.mrao.cam.ac.uk/$\sim$guy/J1118+480/J1118480.list). It is to be noted
that the radio flux ($F_R$) decreases as PCA rate decreases (see, Fig. 3(e) and Fig. 2(a)). So, it appears that the outburst
is totally dominated by the emission from jets/outflows. This motivated us to separate X-ray contribution from the jet/outflow
($F_{ouf}$) from that of the accretion disk or inflowing matter ($F_{inf}$) using the same method as described in JCD17 and Paper-I.

According to TCAF, model normalization parameter is a function of the mass of the black hole ($M_{BH}$), the inclination
angle ($i$) and the distance ($D$) of the system. So when there is no prominent X-ray contribution from other physical
processes (whose effects are not considered in the current version of the TCAF model fits file), for example, jets or outflows
or precession in the disk, $N$ should not vary from one day to another if observed with a given instrument. But during the outburst,
we see a monotonic decrease in $N$ from $10.9$ to $4.41$ (see, Fig. 3(d)). Its variation is roughly similar to $F_R$ (see Fig. 3(e)).
This allows us to conclude that higher $N$ values are possibly required to fit the spectra to match excess contribution of X-rays,
emitted from the base of the jet.

Total X-ray fluxes ($F_X$) are obtained in our PCA spectral analysis band ($3-25$~keV), when all the TCAF model input
parameters (except the mass of the BH) are kept free. In presence of jet or outflow, $F_X$ is the combined contribution of
the radiations from the accretion disk and CENBOL, i.e. from inflowing matter ($F_{inf}$) and from the base of the jet i.e., from
outflowing matter ($F_{ouf}$). Significant X-ray contribution from the jet to the total X-ray influences our spectral fitting
with the TCAF model. Since present version of the TCAF model fits file takes care only the radiation contributed from the
inflowing matter, higher values of normalization are required to fit the spectra as to compensate the excess X-ray radiations
emitted from the jet. In the absence of jet, an almost constant value of model normalization is sufficient to fit the spectra
during the entire outburst (see Molla et al. 2016, 2017; Chatterjee et al. 2016). For the case of 2005 outburst of XTE J1118+480,
we looked into the obtained $N$ values and find that on the last data (MJD=53395.59), the model normalization was at its minimum
value of $4.41$. This implies that on this observation total X-ray flux was contributed only by the radiations from the accretion
disk and CENBOL. So, we may say that on this observation jets X-ray was minimum or negligible. This was indeed observed
via variation of the radio fluxes (see, Fig. 3(e)). So, after confirmation from low $F_R$ at the lowest $N$ observed day,
we may obtain $F_{inf}$ during the entire outburst by refitting spectra with the frozen model $N$ value (at its lowest
observed value, obtained when all model flow parameters are kept free). Now we could separate jet component of the X-ray
fluxes ($F_{ouf}$) in each observations by subtracting $F_{inf}$ from $F_X$ (observed flux when all model flow parameters
are kept free), i.e.,
$$
F_{ouf}=F_X - F_{inf}  \eqno(1)
$$
This method of separating jet contribution of X-ray from total observed X-ray was introduced in Jana et al. (2017).

Although during 2005 outburst of XTE~J1118+480, minimum $N$ value $4.41$ was required to fit spectra, according to Paper-I,
lowest $N$ value was required to fit 2000 outburst spectra was $4.37$. Since Paper-I $N$ value was lowest between two
outbursts, to obtain X-ray contribution only from the disk and CENBOL i.e., inflowing matter, we refitted all 2005 studied
PCA spectra with the TCAF model after keeping model normalization frozen at $4.37$. Finally, $F_{ouf}$ is calculated for each
observation by taking the difference of $F_{inf}$ and $F_X$ as in Eqn (1).
The evolution of $F_X$, $F_{inf}$, $F_{ouf}$, $N$, and $F_R$ are shown in Fig. 3. $F_{ouf}$ shows a similar
variation as of $N$ and $F_R$. When we see the percentage of jet X-ray contribution in total emitted X-ray, as
outburst progresses it decreases (see, Col. 6 of Table 2). Maximum fractional jet X-ray contribution is found
to be $\sim 60\%$.

\begin{figure}[h]
\vskip 0.2cm
\centerline{
           \includegraphics[scale=0.6,width=14.0truecm,angle=0]{fig3.eps}
           }
\caption{Variation of (a) total flux ($F_X$), (b) inflow or disk flux ($F_{inf}$), (c) outflow or jet flux ($F_{ouf}$),
(d) model normalization ($N$), and (e) radio flux ($F_R$) are shown with MJD. $F_X$, $F_{inf}$, $F_{ouf}$ are obtained
TCAF model fitted spectra in $3-25$~keV PCA range in units of $10^{-9}~erg~cm^{-2}~s^{-1}$ unit. $F_R$ of $15~GHz$ Ryle
Telescope is presented in $mJy$ unit and adopted from Brocksopp et al. (2010).}
\label{fig4}
\end{figure}

\section{Summary and Discussions}

The 2005 outburst of XTE J1118+480 (epoch under consideration here) is a shorter duration, lower intense and less studied
compared to its earlier 2000 outburst. To understand the nature of the accretion flow dynamics of the source during this outburst,
we use RXTE PCA and HEXTE combined data for our temporal and spectral analysis. We did not find any prominent signature of
low frequency QPOs in the $0.01$~sec binned Fourier transformed power density spectra, although Shahbaz et al. (2005) reported
a QPO at $\sim 2$~mHz from the optical data in the quiescent state (June 2003), which is in between two outbursts of the source.
The spectral analysis is done using two types of models: the phenomenological power-law (PL) model and physical TCAF solution
based {\it fits} file. Total $21$ observed data of combined PCA and HEXTE spectra in the energy range of $3-100$~keV
(except some observations where due to low s/n, $3-40$~keV data are used) are chosen for spectral fittings from 2005 January 14
(MJD=53384.99) to January 25 (MJD=53395.59). Although triggering of the outburst was reported by Zurita et al. (2005a) on
2005 Jan. 9, RXTE started monitoring it five days later. So we missed an important rapidly evolving rising phase of the outburst.

The evolution of the average PCA rate in $2-25$~keV band implies that on the first observation day (MJD=53384.99), the source
already passed its peak intensity (see, Fig. 2(a)). During the entire period of the outburst, the photon count was found to be
monotonically decreasing. While fitting spectra with the PL model, we noticed the same behaviour for the PL flux (see Fig. 2(b)).
We tried to fit the spectra with combined DBB plus PL models, but `ftest' suggest that DBB component is insignificant. This
means that during the entire outburst, non-thermal PL photons are highly dominating over the low or insignificant thermal disk
component. PL model fitted photon index ($\Gamma$) values are obtained in a lower range ($\sim1.75-1.87$), implies that during
the entire phase of the outburst, the source was in the hard state. If we relook at the $\Gamma$ values, we see an increasing
trend of it (from $1.81$ to $1.87$) i.e., rise in softness in the late declining phase of the outburst. There is no physical
explanation of this $\Gamma$ from the PL model fitted spectral analysis.

During the entire outburst, high dominance of the sub-Keplerian halo rate ($\dot m_h$) over the Keplerian disk rate ($\dot m_d$)
in the presence of strong shock ($R>3$) far away from the BH ($X_s>245~r_s$) are observed. But if we look into the evolution
of accretion rates, in the late declining phase (specifically in last observations) a slow rise in $\dot m_d$ is observed although
$\dot m_h$ is decreased. This explains the gradual softening i.e., increase in $\Gamma$ in the late declining phase of the outburst.
This feature is quite uncommon in an outburst of a transient BHC. This was due to the vanishing of the jet effect.
A similar late softening was also observed due to the slow rise in $\dot m_d$ during the late declining phase of 2000 outburst
(Paper-I). So the source has been consistently behaving in a non-conventional way even after a gap of five long years.
This points to the fundamental configuration of the binary and compactness of the system.

Radio flares are believed to be a tracer of jets or outflows. Similar to the 2000 outburst, high radio fluxes are also observed
during the current outburst of XTE~J1118+480. $F_R$ showed similar declining variation as of PCA rate, PL flux, $F_X$.
Now, when we looked at the TCAF model fitted $N$ values, similar nature of monotonically decreasing values (from $10.9$
to $4.41$) are observed. This means that higher $N$ were required to fit spectra when stronger jets were present.
We estimated $F_{inf}$ by refitting the spectra with $N$ frozen at its minimum observed value ($=4.37$) during its 2000 and 2005
outbursts (when spectral fits are done with all TCAF model parameters as free).
We estimate the X-ray flux contribution from jets ($F_{ouf}$) and find that the fractional jet X-ray contribution ($F_{ouf}$/$F_X$)
is maximum $\sim 60\%$. Evolution of $F_{ouf}$ looks similar to $F_R$. It is to note that the outburst is suppressed with the
decrease in $F_{ouf}$. Therefore, similar to the 2000 outburst of XTE~J118+480, this so-called `failed outburst' could be termed as
jet activity dominated outburst.

\section*{Acknowledgements}
This work made use of PCA and HEXTE data of NASA's RXTE satellite.
Research of D.D. and S.K.C. is supported in part by the Higher Education Dept. of the Govt. of West Bengal, India.
D.D. and S.K.C. also acknowledge partial support from ISRO sponsored RESPOND project (ISRO/RES/2/418/17-18) fund.
D.C. and D.D. acknowledge support from DST/SERB sponsored Extra Mural Research project (EMR/2016/003918) fund.
A.J. and D.D. acknowledge support from DST/GITA sponsored India-Taiwan collaborative project (GITA/DST/TWN/P-76/2017) fund.
A.J.  acknowledges CSIR SRF fellowship (09/904(0012)2K18 EMR-1).
K.C. acknowledges support from DST/INSPIRE fellowship (IF170233).

\bibliographystyle{raa}