Effective implementation of $l_0$-Regularised Compressed Sensing with Chaotic-Amplitude-Controlled Coherent Ising Machines

TL;DR

This paper demonstrates that a closed-loop chaotic-amplitude-controlled coherent Ising machine improves the accuracy and effectiveness of $l_0$-regularised compressed sensing compared to open-loop systems, using optical and MRI data.

Contribution

It introduces a closed-loop CIM with chaotic amplitude control for enhanced $l_0$-regularised compressed sensing, surpassing open-loop system performance.

Findings

Improved accuracy over open-loop systems.

Wider effectiveness range demonstrated.

Validated with artificial and MRI data.

Abstract

Coherent Ising Machine (CIM) is a network of optical parametric oscillators that can solve large-scale combinatorial optimisation problems by finding the ground state of an Ising Hamiltonian. As a practical application of CIM, Aonishi et al., proposed a quantum-classical hybrid system to solve optimisation problems of $l_0$-regularisation-based compressed sensing. In the hybrid system, the CIM was an open-loop system without an amplitude control feedback loop. In this case, the hybrid system is enhanced by using a closed-loop CIM to achieve chaotic behaviour around the target amplitude, which would enable escaping from local minima in the energy landscape. Both artificial and magnetic resonance image data were used for the testing of our proposed closed-loop system. Compared with the open-loop system, the results of this study demonstrate an improved degree of accuracy and a wider range of effectiveness.

Equations

$$\hat{x}=\underset{x\in {\mathbb{R}}^N}{\mathrm{argmin}}\Vert x{\Vert }_{p}subjecttoy=Ax.$$

$$(\hat{R},\hat{\sigma })=\underset{\sigma \in \{0,1{\}}^N}{\mathrm{argmin}}\underset{R\in {\mathbb{R}}^N}{\mathrm{argmin}}\left(\Vert y-A(\sigma \circ R){\Vert }_2^2\right)subjectto\Vert \sigma {\Vert }_0\le \Omega .$$

$$(R,\sigma )=\underset{\sigma \in \{0,1{\}}^N}{\mathrm{argmin}}\underset{R\in {\mathbb{R}}^N}{\mathrm{argmin}}\left(\frac{1}{2}\Vert y-A(\sigma \circ R){\Vert }_2^2+\lambda \Vert \sigma {\Vert }_0\right).$$

$$\text{mathbcal}H=\sum _{r<r^{\prime }}^N\sum _{k=1}^M{A}_r^k{A}_{r^{\prime }}^k{R}_r{R}_{r^{\prime }}{\sigma }_r{\sigma }_{r^{\prime }}-\sum _{r=1}^N\sum _{k=1}^M{y}^k{A}_r^k{R}_r{\sigma }_r+\lambda \sum _{r=1}^N{\sigma }_r,$$

$${\left(\frac{d{c}_r}{dt}\right)}_{inj,r}=\left(\text{absolutevalue}{h}_r-\eta \right).$$

$$\text{centering}{h}_r=-\sum _{r^{\prime }=1(\ne r)}^N\sum _{k=1}^M{A}_r^k{A}_{r^{\prime }}^k{R}_{r^{\prime }}H(c_{r^{\prime }})+\sum _{k=1}^M{A}_r^k{y}^k,\text{@add@centering}$$

$$h_r=-\sum _{r^{\prime }=1(\ne r)}^N\sum _{k=1}^M{A}_r^k{A}_{r^{\prime }}^k{R}_{r^{\prime }}H\left(c_{r^{\prime }}\right)+\sum _{r^{\prime }=1}^N\sum _{k=1}^M{A}_r^k{A}_{r^{\prime }}^k{x}_{r^{\prime }}{\xi }_{r^{\prime }},$$

$${\left(\frac{d{\mu }_r}{dt}\right)}_{inj,r}=j{e}_r\left(R_r{h}_r-\frac{{\eta }^2}{4}\sqrt{\frac{\tau }{g^2}}\right),$$

$$\text{centering}\frac{d}{dt}{e}_r=-\beta \left(g^2{\tilde{\mu }}_r^2-\tau \right)e_r,\text{@add@centering}$$

$${\tilde{\mu }}_r={\mu }_r+\sqrt{\frac{1}{4j}}{W}_{R,r},$$

$$\text{centering}{h}_r=-\sum _{r^{\prime }=1(\ne r)}^N\sum _{k=1}^M{A}_r^k{A}_{r^{\prime }}^k{R}_{r^{\prime }}\frac{1}{2}\left({\tilde{\mu }}_{r^{\prime }}+\sqrt{\frac{\tau }{g^2}}\right)+\sum _{k=1}^M\sqrt{\frac{\tau }{g^2}}{A}_r^k{y}^k,\text{@add@centering}$$

$$h_r=-\sum _{r^{\prime }=1(\ne r)}^N\sum _{k=1}^M{A}_r^k{A}_{r^{\prime }}^k{R}_{r^{\prime }}\frac{1}{2}\left({\tilde{\mu }}_{r^{\prime }}+\sqrt{\frac{\tau }{g^2}}\right)+\sum _{r^{\prime }=1}^N\sum _{k=1}^M\sqrt{\frac{\tau }{g^2}}{A}_r^k{A}_{r^{\prime }}^k{x}_{r^{\prime }}{\xi }_{r^{\prime }}.$$

$$\begin{multlined}\frac{d}{dt}c_{r}=\left[-1+p-{\left(c_{r}^{2}+s_{r}^{2}\right)}\right]c_{r}+\widetilde{K}\left(\dfrac{dc_{r}}{dt}\right)_{inj,r}+\\ {g}\sqrt{\left(c_{r}^{2}+s_{r}^{2}\right)+\frac{1}{2}}W_{r,1},\end{multlined}\frac{d}{dt}c_{r}=\left[-1+p-{\left(c_{r}^{2}+s_{r}^{2}\right)}\right]c_{r}+\widetilde{K}\left(\dfrac{dc_{r}}{dt}\right)_{inj,r}+\\ {g}\sqrt{\left(c_{r}^{2}+s_{r}^{2}\right)+\frac{1}{2}}W_{r,1},$$

$$\begin{multlined}\frac{d}{dt}s_{r}=\left[-1-p-{\left(c_{r}^{2}+s_{r}^{2}\right)}\right]s_{r}+{g}\sqrt{\left(c_{r}^{2}+s_{r}^{2}\right)+\frac{1}{2}}W_{r,2}.\end{multlined}\frac{d}{dt}s_{r}=\left[-1-p-{\left(c_{r}^{2}+s_{r}^{2}\right)}\right]s_{r}+{g}\sqrt{\left(c_{r}^{2}+s_{r}^{2}\right)+\frac{1}{2}}W_{r,2}.$$

$$\begin{multlined}\dfrac{d}{dt}\mu_{r}=-\left(1-p+j\right)\mu_{r}-g^{2}\mu_{r}^{3}+\sqrt{j}\left(V_{r}-\frac{1}{2}\right)W_{R,r}+{{K}}\left(\frac{d\mu_{r}}{dt}\right)_{inj,r},\end{multlined}\dfrac{d}{dt}\mu_{r}=-\left(1-p+j\right)\mu_{r}-g^{2}\mu_{r}^{3}+\sqrt{j}\left(V_{r}-\frac{1}{2}\right)W_{R,r}+{{K}}\left(\frac{d\mu_{r}}{dt}\right)_{inj,r},$$

$$\begin{multlined}\dfrac{d}{dt}V_{r}=-2\left(1-p+j\right)V_{r}-6g^{2}\mu_{r}^{2}V_{r}+1+j+2g^{2}\mu_{r}^{2}-2j\left(V_{r}-\frac{1}{2}\right)^{2}.\end{multlined}\dfrac{d}{dt}V_{r}=-2\left(1-p+j\right)V_{r}-6g^{2}\mu_{r}^{2}V_{r}+1+j+2g^{2}\mu_{r}^{2}-2j\left(V_{r}-\frac{1}{2}\right)^{2}.$$

$$\begin{multlined}\dfrac{d}{dt}\mu_{r}=-\left(1-p+j\right)\mu_{r}-{g^{2}\mu_{r}\left(\mu_{r}^{2}+2n_{r}+m_{r}\right)}+\sqrt{j}\left(m_{r}+n_{r}\right)W_{R,r}\\ +{{K}}\left(\frac{d\mu_{r}}{dt}\right)_{inj,r},\end{multlined}\dfrac{d}{dt}\mu_{r}=-\left(1-p+j\right)\mu_{r}-{g^{2}\mu_{r}\left(\mu_{r}^{2}+2n_{r}+m_{r}\right)}+\sqrt{j}\left(m_{r}+n_{r}\right)W_{R,r}\\ +{{K}}\left(\frac{d\mu_{r}}{dt}\right)_{inj,r},$$

$$\begin{multlined}{\dfrac{d}{dt}n_{r}=-2\left(1+j\right)n_{r}+2pm_{r}-2g^{2}\mu_{r}^{2}\left(2n_{r}+m_{r}\right)}{-j\left(m_{r}+n_{r}\right)^{2}},\end{multlined}{\dfrac{d}{dt}n_{r}=-2\left(1+j\right)n_{r}+2pm_{r}-2g^{2}\mu_{r}^{2}\left(2n_{r}+m_{r}\right)}{-j\left(m_{r}+n_{r}\right)^{2}},$$

$$\begin{multlined}{\dfrac{d}{dt}m_{r}=-2\left(1+j\right)m_{r}+2pn_{r}-2g^{2}\mu_{r}^{2}\left(2m_{r}+n_{r}\right)+}p\\ -g^{2}\left(\mu_{r}^{2}+m_{r}\right)-j\left(m_{r}+n_{r}\right)^{2}.\end{multlined}{\dfrac{d}{dt}m_{r}=-2\left(1+j\right)m_{r}+2pn_{r}-2g^{2}\mu_{r}^{2}\left(2m_{r}+n_{r}\right)+}p\\ -g^{2}\left(\mu_{r}^{2}+m_{r}\right)-j\left(m_{r}+n_{r}\right)^{2}.$$

$${\sigma }_r=Heaviside\left(x_r\right)=\{\begin{array}{ll}1,& \left(x_r>0\right)\\ 0,& \left(x_r\le 0\right).\end{array}$$

$$R_r\sum _{k=1}^M{\left(A_r^k\right)}^2={\sigma }_r{\mathbb{H}}_r,$$

$${\mathbb{H}}_r=-\sum _{r^{\prime }=1(\ne r)}^N\sum _{k=1}^M{A}_r^k{A}_{r^{\prime }}^k{R}_{r^{\prime }}{\sigma }_{r^{\prime }}+\sum _{k=1}^M{A}_r^k{y}^k.$$

$$p=(p_{thr}-d)+\frac{2d}{1+e^{-\left(\frac{t-4}{2}\right)}}.$$

$$p=1.5\times {\left(\frac{t}{5}\right)}^2.$$

$${\eta }_i=\max \left[{\eta }_{init}\left(1-\frac{i}{velo}\right),{\eta }_{end}\right].$$

$$\left[\begin{array}{c}{y}^1\\ y^2\\ ⋮\\ y^M\end{array}\right]=\left[\begin{array}{cccc}{A}_1^1& A_2^1& \cdots & A_N^1\\ A_1^2& A_2^2& \cdots & A_N^2\\ ⋮& ⋮& \ddots & ⋮\\ A_1^M& A_2^M& \cdots & A_N^M\end{array}\right]\left[\begin{array}{c}{\xi }_1{x}_1\\ {\xi }_2{x}_2\\ ⋮\\ {\xi }_N{x}_N\end{array}\right]+\left[\begin{array}{c}{w}_{noise}^1\\ w_{noise}^2\\ ⋮\\ w_{noise}^M\end{array}\right].$$

$$x=\mathrm{argmin}(\Vert y-SFx{\Vert }_2^2+\frac{1}{2}\gamma \Vert {\Delta }_{v}x{\Vert }_2^2+\frac{1}{2}\gamma \Vert {\Delta }_{h}x{\Vert }_2^2+\lambda \Vert \Psi x{\Vert }_0).$$

$$h^z=SF{\Psi }^{T}y,$$

$$\tilde{J}=\Psi F^T{S}^{T}SF{\Psi }^T+\gamma \Psi {\Delta }_v^T{\Delta }_v{\Psi }^T+\gamma \Psi {\Delta }_h^T{\Delta }_h{\Psi }^T.$$

Full text

1 Introduction

Compressed sensing (CS) is a method of reconstructing a high-dimensional signal or image based on highly downsampled measurements.

There has been considerable interest in it across a wide range of fields and applications. Such as in the field of astronomy, a possible way to transmit data to Earth from spacecraft CSastro has been attempted. And there are proposed methods with CS on astronomical image compression and in compression on remotely sensed data CSastro2 ; CSastro3 ; CSastro4 as well. And in radar technologies for the reconstruction of the target image CS has been used CSsig2 . On the other hand in the medical field using embedded compression using CS to improve energy efficiency in Electrocardiogram (ECG) machines has been proposed CSbiosig .

[TABLE]

The above equation shows an observed signal $y\in\mathbb{R}^{M}$, an observation matrix $A\in\mathbb{R}^{M\times N}$, and a source signal $x\in\mathbb{R}^{N}$. Hereafter, the ratio of the number of non-zero entries in $x$ to $N$ is defined as the sparseness $a$, and the ratio of $M$ to $N$ is defined as the compression ratio $\alpha$. Since $l_{1}$-norm CS is a convex optimisation problem, there are many efficient algorithms for optimisation of $l_{1}$-norm CS that are widely applied in the real-world problems mentioned above. However, there has been a suggestion that $l_{0}$-norm CS should outperform $l_{1}$-norm CS since the $l_{1}$-norm penalty does not lead to any solution shrinkage obuchi ; kabashima . In the thermodynamic limit $N$, $M$ $\longrightarrow$ $\infty$ with $\alpha=M/N$ kept fixed, an $l_{0}$-norm CS’s threshold for $a$, determining whether or not the problem has a solution with no error, is larger than that of $l_{1}$-norm CS’s obuchi ; kabashima . Nonetheless, the optimisation in $l_{0}$-norm CS is challenging since it involves combinatorial optimisation.

Numerous attempts have been made to overcome the issue in $l_{0}$-norm CS optimisations. $l_{0}$-norm CS can be formulated as a two-fold optimisation twofold1 ; twofold2 .

[TABLE]

Here $R\in\mathbb{R}^{N}$ and $\sigma\in\left\{{0,1}\right\}^{N}$ correspond to the source signal and support vector, respectively. Especially, each entry in the support vector taking either 0 or 1 represents whether each entry in the source signal is zero or non-zero. The condition $\|\sigma\|_{0}\leq\Omega$ is a sparsity-inducing prior for constraining the number of non-zero entries to be $\Omega$. Therefore, the optimisation with respect to $\sigma$ can be regarded as a quadratic-constrained binary optimisation problem to find a ground state of a two-state Potts Hamiltonian. Based on this formulation, simulated annealing (SA) algorithm has been attempted obuchi . On the other hand, Aonishi et al., attempted to solve optimisation problems of $l_{0}$-norm CS with a quantum-classical hybrid approach. $l_{0}$-norm CS implemented with the hybrid system is given as a regularisation form as follows Aonishi .

[TABLE]

The element-wise representation of Eq. (3) gives the following Hamiltonian.

[TABLE]

where an element $A^{k}$ in $A$, an element $y^{k}$ in $y$, an element $R_{r}$ in $R$ and an element $\sigma_{r}$ in $\sigma$. Optimisation with respect to $\sigma$ in Eq. (4) is a quadratic unconstrained binary optimisation (QUBO) problem, which is implementable with a quantum machine such as the coherent Ising machine (CIM) Aonishi ; cimqubo1 ; cimqubo2 ; cimqubo3 . In the quantum-classical hybrid approach to conducting $l_{0}$-regularised CS, $\sigma$ is optimised by the CIM while $R$ is optimised by a Classical Digital Processor (CDP) (see Fig. 1).

The CIM architecture in the hybrid approach was an open-loop (OL) CIM with the Zeeman term. The hybrid approach with the OL-CIM is hereafter referred to as OL-CIM-CDP. Note that the OL means the lack of feedback loop for amplitude control described below. It has been reported that the imbalance in the size of the interaction term and the Zeeman term degrades the system performance AonishiCDMA . To balance these terms, for the local field, the measured-amplitudes were binarised. OL-CIM-CDP in this formulation outperformed SA on the regularisation form Aonishi .

The close-loop CIM, in which the amplitudes of optical parametric oscillator (OPO) pulses are controlled to a target value, have been proposed to improve the performance of CIM’s ground-state search kako ; Inui2022 . Especially, introducing auxiliary nonlinear dynamics forcefully trying to equalise to a target value results in chaotic behaviour around the target in the CIM which may result in escaping from local minima in the energy landscape. This chaotic method is referred to as chaotic amplitude control (CAC) AmpLeleu ; Leleu1 ; sam ; kako ; Inui2022 . Recently, Inui et al., have proposed an approach to efficiently incorporate the Zeeman terms in CAC-CIM by scaling the Zeeman terms with target amplitude to match that of the interaction term Inui2022 .

In this paper, following Inui et al.’s approach, we modify the CAC-CIM for performing QUBO in $l_{0}$-regularised CS and attempt to improve the performance of the hybrid CIM-CDP system by replacing the OL-CIM with the CAC-CIM with the Zeeman term (see Fig. 1). The hybrid system proposed here is hereafter referred to as CAC-CIM-CDP. Firstly, to demonstrate the effectiveness of CAC-CIM for performing QUBO in the support estimation, we compare the performance of CAC-CIM to those of OL-CIM and SA. Then, to demonstrate the effectiveness of CAC-CIM-CDP for performing an alternating minimisation, we compare the performance of CAC-CIM-CDP to that of OL-CIM-CDP on artificial random data, as well as magnetic resonance imaging (MRI) data.

2 Results

2.1 Alternating minimisation algorithm

Alternating minimisation procedures on CAC-CIM-CDP and OL-CIM-CDP are summarised in Algorithm 1 and Algorithm 2, respectively. This type of minimisation suggests the back-and-forth optimisation performed between the CIM and CDP. CIM passes the optimisation results to the CDP after optimising the support, as shown in Fig. 1. The CDP then optimises the signal and sends the resulting signal to the CIM for support optimisation. In Algorithm 1 and Algorithm 2, indicate the number of iterations of alternating minimisation, the initial values and the integration interval for stochastic differential equations (SDEs) of CIM and so on. The schedules of the pump rate, threshold and target amplitude are given in section 4.3.

2.2 Outline of the CIM models and injection field for QUBO on support estimation

On CIM, $l_{0}$-regularised CS is performed by updating the injection field dictated by the local field, which is determined by the gradient of the QUBO Hamiltonian Eq. (4) with respect to the spin coordinates. Aonishi et al., proposed OL-CIM-CDP, which is based on an open-loop injection scheme Aonishi . They used the CIM model expressed as the Wigner stochastic differential equation (W-SDE) Eq. (13) and Eq. (14) (in Methods) with the following injection field.

[TABLE]

[TABLE]

Here, $h_{r}$ is the local field expressed as Eq. (6). $R_{r}$ is the signal value estimated by the CDP. $c_{r}$ is the in-phase amplitude of the $r$-th OPO pulse, and $H(c_{r})$ is the binarised in-phase amplitude by the Heaviside step function as proposed in the discrete simulated bifurcation disSimBif . $\eta$ is the threshold which is related to the $l_{0}$-regularisation parameter $\lambda$ by $\eta=\sqrt{2\lambda}$ according to the Maxwell rule (see Aonishi for a detailed explanation). In the local field Eq. (6), the mutual interaction is $\tilde{J}_{rr^{\prime}}=-\sum_{k=1}^{M}A_{r}^{k}A_{r^{\prime}}^{k}$ and the Zeeman term is $\sum_{k=1}^{M}A_{r}^{k}y^{k}$. Substituting the observation model Eq. (26) (in Section 4.4) into Eq. (6) when $w_{noise}=0$ (no observation noise), the local field Eq. (6) can be expressed as follows.

[TABLE]

where $x_{r}$ is the true signal value, $\xi_{r}$ is the true support taking 1 or 0. The Zeeman term in the second term of Eq. (7) can be regarded as the matched filter, in which $A^{T}A$ is calculated. The mutual interaction term in the first term plays a role in removing off-diagonal elements ($r\neq r^{\prime}$) corresponding to cross-talk noise in the Zeeman term, which are induced by the cross-correlation among the column vectors $A_{1},...,A_{N}$ in $A$. To obliterate the cross-talk noise, the in-phase amplitude $c_{r}$ needs to be the same as the amplitude of $\xi_{r}$ if $R_{r}=x_{r}$. Hence, $c_{r}$ is binarised to either 1 or 0. In Fig. 2e, a typical evolution of $c_{r}$ in the open-loop-type W-SDE is illustrated. $c_{r}$ does not keep the same amplitude as that of $\xi_{r}$ and increases with increasing the pump rate.

In this paper, we propose CAC-CIM-CDP, based on a closed-loop injection scheme with CAC. The idea of CAC for CIM was first introduced by Leleu et al., AmpLeleu . It simply states that forcefully trying to equalise the amplitudes of the system to a specific value (in CAC, target amplitude $\tau$) may result in a chaotic behaviour in the system which may result in escaping from local minima in the energy landscape. In this paper, we used two CIM models expressed as W-SDE Eq. (15) and (16) and Positive-$P$ stochastic differential equation (P-SDE) Eq. (17)-(19) (in Section 4.1) commonly having the following injection field with CAC feedback.

[TABLE]

[TABLE]

[TABLE]

[TABLE]

where $h_{r}$ is the local field expressed as Eq. (11), $e_{r}$ is the auxiliary variable for the error feedback in the CAC feedback loop, and $\tau$ indicates the target amplitude for the CAC. $R_{r}$ is the signal value estimated by the CDP, which is the same as that of OL-CIM-CDP. $\eta$ is the threshold given by $\eta=\sqrt{2\lambda}$, which is introduced to keep consistency with OL-CIM-CDP. As described in Section 4.1 in Methods, $j$ is the normalised out-coupling rate for optical homodyne measurement, and $g^{2}$ is the nonlinear saturation parameter of the CIM which determines the abrupt jump of the photon number at the OPO threshold and the amplitude of the quantum noise present in CIM. $\tilde{\mu}_{r}$ implies the measured-amplitude, and $w(R,r)$ is the independent real Gaussian noise process, which is the same as that in W-SDE (15) and P-SDE (17). In the local field Eq. (11), the mutual interaction is $\tilde{J}_{rr^{\prime}}=-\sum_{k=1}^{M}A_{r}^{k}A_{r^{\prime}}^{k}$ and the Zeeman term is $h_{r}^{z}=\sqrt{{\tau/g^{2}}}\sum_{k=1}^{M}A_{r}^{k}y^{k}$. Substituting the observation model Eq. (26) into Eq. (11) when $w_{noise}=0$ (no observation noise), the local field Eq. (11) can be expressed as follows.

[TABLE]

In Fig. 2a and 2b, the typical evolution of normalised measured-amplitude $g\tilde{\mu}_{r}$ are shown. The corresponding error evolution is indicated in Fig. 2c and 2d. Due to the CAC feedback loop, as shown in Fig. 2a and 2b, if the squared-amplitude of DOPO is smaller than $\tau$, $e_{r}$ exponentially increases and vice-versa, and the measured-amplitude $\tilde{\mu}_{r^{\prime}}$ is maintained around $\sqrt{{\tau/g^{2}}}$. Therefore, because $1/2(\tilde{\mu}_{r^{\prime}}+\sqrt{{\tau/g^{2}}})$ in Eq. (12) can take around [math] or $\sqrt{{\tau/g^{2}}}$, the mutual interaction term and the Zeeman term scales are balanced, and crosstalk noise, i.e. off-diagonal elements, is eliminated from the Zeeman term as described in OL-CIM-CDP. Moreover, as shown in Fig. 2a and 2b, it is important to note that intermediate solutions are destabilised. By doing so, CAC introduced CIM is able to keep searching for an answer until the maximum run-time has been reached. By taking the support vector that is generated by CIM at the end of each trajectory, we are evaluating the solution to estimate the support for the simulations in this paper.

2.3 Comparison with Simulated Annealing

Here our purpose is to demonstrate that CAC feedback is effective on CIM by comparing CAC-CIM to OL-CIM and SA. We follow the Metropolis algorithm for $l_{0}$-regularised CS stated in Aonishi . As same as in Aonishi , 1000 samples of the observation matrix and source signal and true support vector are randomly generated according to Section 4.5 under $N=500$, $\alpha=a=0.6$, $w_{noise}=0$ (no observation noise). With the same observation matrices, source signals, and support vectors in all models, we statistically evaluate how well CAC-CIM estimates support in comparison to OL-CIM and SA when all $R_{r}$ are fixed to be the source signal $x_{r}$. To measure the support estimation quality, we used the direction cosine defined as ${\sum_{r=1}^{N}\xi_{r}\sigma_{r}}/{\sqrt{\sum_{r=1}^{N}\xi_{r}\sum_{r=1}^{N}\sigma_{r}}}$ where $\left(\xi_{1},...,\xi_{N}\right)$ is the true support vector and $\left(\sigma_{1},...,\sigma_{N}\right)$ is the estimated one. When the estimation is perfect, the direction cosine is equal to 1. We selected $\eta=0.05$ corresponding to $l_{0}$-regularisation parameter $\lambda=\eta^{2}/2=0.00125$ as in Aonishi .

First, we evaluate the temporal profiles of the optimisation processes for the support estimation in CAC-CIM (Wigner), CAC-CIM (Positive-$P$), OL-CIM and SA. The upper three graphs (from left to right, CAC-CIM (Wigner), CAC-CIM (Positive-$P$) and OL-CIM respectively) in Fig. 3a show the change in the direction cosine of the three CIM models depending on the runtime on the CPU and the wall-clock time of physical CIM. For CAC-CIM (Wigner), and CAC-CIM (Positive-$P$) models, $20\times$ photon’s lifetimes of integral interval (with 1000 time-steps) for the SDEs are about 105ms and 68ms of run-time respectively, and for OL-CIM, $5\times$ photon’s lifetime of integral interval (with 50 time-steps) for the SDE is about 11ms. The physical CIM’s wall-clock time for this optimisation is roughly estimated to be around 0.5ms, which can be estimated from the round-trip time of $N=500$ and the time-steps-to-solution for the Sherrington-Kirkpatrick problem with $N=500$ sam . The direction cosine of these CIM models converged to about 1 by these run-times. The lower two graphs in Fig. 3a show the change in the direction cosine of SA depending on the runtime on CPU under constant temperature at $T=0$ and exponential cooling scheduling from $T=0.02$ to $0.00002$. We adjusted the Monte-Carlo steps of SA (bottom two graphs of Fig. 3a) to accompany the wall-clock time of physical CIM (0.5ms) and the run-time of CAC-CIM (Wigner) (105ms). In our computational environment, the number of Monte Carlo steps for SA with runtimes of 0.5ms and 105ms is about 230 and 46000 steps, respectively. In SA, the direction cosine converged to about 1 by 105ms, while that did not by 0.5ms.

Next, we compare the histogram of the final states of direction cosines in CAC-CIM (Wigner), CAC-CIM (Positive-$P$), OL-CIM and SA. The upper three graphs in Fig. 3b indicate the histogram of the three CIM models (CAC-CIM (Wigner), CAC-CIM (Positive-$P$), OL-CIM, respectively), while the lower two graphs in Fig. 3b show the histograms of SA at run-times of 0.5ms and 105ms under zero temperature and exponential cooling schedules respectively. Comparing these graphs, the proportion of the direction cosines of CAC-CIM (Wigner) and CAC-CIM (Positive-$P$) close to 1 is higher than those of OL-CIM and SA. The two-sample one-sided Kolmogorov-Smirnov test suggests that the histograms of the final direction cosines of CAC-CIM (Wigner) and CAC-CIM (Positive-$P$) are significantly biased toward 1 compared with all of those of OL-CIM and SA (P-value $<$ 0.0001).

The above results thus demonstrate that CAC-CIM outperformed OL-CIM on support vector estimation and outperformed SA within the same run-time.

2.4 Comparison with ground state predicted with statistical mechanics on alternating minimisation

We compare CAC-CIM-CDP’s capability to find the ground state with that of OL-CIM-CDP. In our previous study, we derived the macroscopic parameter equation (Eq. (26)-(28) in Aonishi ) using a non-equilibrium statistical mechanics method to show the performance limit of OL-CIM-CDP. In the limit of the saturation parameter $g^{2}\rightarrow 0$, the macroscopic parameter equation derived in the previous study is consistent with that for a two-state Potts spin system defined by the QUBO Hamiltonian Eq. (4). Therefore, the macroscopic parameter equation in this limit can predict the ground state of the Hamiltonian. Through a comparison of solutions of CAC-CIM-CDP and OL-CIM-CDP with a solution of the macroscopic parameter equation in the limit of $g\rightarrow 0$, we demonstrate the efficacy of CAC feedback on the alternating minimisation for optimising the Hamiltonian.

The precondition for applying statistical mechanics is that the values of all entries in the observation model Eq. (26), which is the premise of Eq. (3) and Eq. (4), are randomly determined as described in Section 4.5. To compare solutions of the models with the ground state predicted with statistical mechanics, 10 samples of the observation matrix and source signal and true support vector are randomly generated according to Section 4.5 under $N=2000$ and various values of $a,\alpha$ and $\nu$. Here $\nu$ indicates the standard deviation of the observation noise ($w_{noise}$). Then, we execute Algorithms 1 and 2 for the alternating minimisation in CAC-CIM-CDPs (Wigner and Positive-$P$) and OL-CIM-CDP sharing the same samples of observation matrices, source signals and support vectors. Here for Fig. 4, $\eta_{init}=0.6$ and $\eta_{init}=0.8$ was used for CAC-CIM-CDP models and OL-CIM-CDP respectively. $\eta_{end}$ was set to 0.18 in Fig. 4a and Fig. 4b while in Fig. 4c and Fig. 4d $\eta_{end}$ was set to 0.35.

The marks in Fig. 4 show the averaged root-mean-square-error (RMSE) calculated as $\sqrt{1/N\sum_{r=1}^{N}\left(R_{r}\sigma_{r}-x_{r}\xi_{r}\right)^{2}}$ of sampled solutions obtained from OL-CIM-CDP, Wigner and Positive-$P$ of CAC-CIM-CDPs. Here $\sigma_{r}$ is calculated as stated in Eq. (20). The black solid lines in Fig. 4 indicate RMSE at the ground state corresponding to successful signal retrieval, which is predicted with statistical mechanics. RMSEs of Wigner and Positive-$P$ CAC-CIM-CDPs tend to keep a better consistency with that of the ground state compared to OL-CIM-CDP for various values of $a,\alpha$ and $v$. Especially as shown in Figs. 4b and 4d, RMSE of OL-CIM-CDP tend to deviate gradually from that of the ground state as increasing $a$, while both Wigner and Positive-$P$ CAC-CIM-CDPs keep up a better consistency with the theoretical prediction.

2.5 Application to Sparse MRI

We evaluate the performance of CAC-CIM-CDP, OL-CIM-CDP and LASSO Tibshirani on MRI data.

In the following numerical experiment, we used two different-sized sparse images ($64\times 64$ and $128\times 128$ pixels) spanned by a Haar basis function. Detailed explanations of the two images we used as the source images are given in Section 4.6 in Methods. In accordance with our previous work Aonishi , we sought to reconstruct the two images from the undersampled $k$-space data and by solving the optimisation problem defined in Eq. (27) (see Section 4.6). To realise the optimisation problem in Eq. (27) on CIM, the Haar wavelet transform coefficients are estimated with the mutual interaction term and the Zeeman term constructed according to Eq. (28) and (29) in Section 4.6. The compression rate of the $k$-space data from the $64\times 64$ and $128\times 128$ images is 0.4 and 0.3 respectively. And the sparseness of the images is 0.212 and 0.178 respectively. As the solver for CDP, we used the Conjugate Gradient Descent method (further details on CDP optimisation refer to Section 4.2).

In Fig. 5a and Fig. 5b, for 10 simulations the average RMSE value is indicated for each threshold $\eta$ for $64\times 64$ and $128\times 128$ images respectively. As for the minimum RMSE in the $64\times 64$ case, LASSO (black line), OL-CIM-CDP (red), CAC-CIM-CDP (Wigner) (green) and CAC-CIM-CDP (Positive-$P$)’s (blue) can be stated as, 0.0292, 0.0216, 0.0182 and 0.0182 respectively (for the corresponding reconstructions see Fig. 6). In the $128\times 128$ case, the minimum RMSE is 0.0276, 0.0242, 0.0209 and 0.0209 respectively (for the corresponding reconstructions see Fig. 7). Comparing the RMSE values acquired it is clear that CAC-CIM-CDP models have a better average performance compared to the other approaches in both image sizes. And even after reaching the optimal reconstruction for the given parameters, CAC-CIM-CDP tends to keep up a minimal error rate compared to LASSO and OL-CIM-CDP. This indicates that the effective range of CAC-CIM-CDP is much wider than OL-CIM-CDP. In both image sizes, the Wigner and Positive-$P$ variations of CAC-CIM-CDP produce identical RMSE results.

In Fig. 6 and 7 the minimal RMSE constructions are shown for LASSO, OL-CIM-CDP, CAC-CIM-CDP (Wigner) and CAC-CIM-CDP (Positive-$P$). In Fig. 7, only CAC-CIM-CDP (Positive-$P$)’s reconstruction is shown because it is clear that both CAC-CIM-CDP (Wigner) and CAC-CIM-CDP (Positive-$P$)’s performance is identical. In the $64\times 64$ image reconstruction when RMSE values are compared, CAC-CIM-CDP models have better reconstruction accuracy. The enlarged portions indicate the difference in pixel identification of each model compared to the initial resized image. Considering both simulations it is clear that even though the system size increases, proposing models have the upper hand in performing an accurate reconstruction compared to other models.

3 Discussion

In this paper, we have proposed an improved CIM approach to solve $l_{0}$-regularised compressed sensing problems. The proposed algorithm has shown that it can outperform the previously proposed algorithm accuracy-wise in all the simulations performed. With the OL-CIM algorithm, the CIM model in use was lacking the CAC feedback for chaotically exploring solutions. Therefore, CAC-CIM has been able to provide convergence to a better solution than OL-CIM. One factor to emphasise here is that CAC does not guarantee convergence to the ground state. Even the ground state is reached, due to the forceful equalisation to $\tau$ may prevent from stopping there. Even though this is the case in this paper, CAC has been shown to be effective especially when the problem instances are relatively harder in both artificial random data and MRI data.

3.1 Effect of system size on performance

The introduction of CAC has previously been shown to have better performance with small-scale frustrated Ising problem instances Inui2022 . In this manuscript, we have demonstrated the applicability of CAC for real-world combinatorial optimisation problems (in this case Compressed sensing) where the problem instances with a Zeeman term are mapped to a QUBO formulation that is large-scale. The simulations with random artificial data on various system sizes are illustrated in Supplementary note 1. Even though the performance increase is present, in very large system sizes such as in $128\times 128$, it is clear that the RMSE gap between CAC-CIM-CDP and OL-CIM-CDP is smaller compared to $64\times 64$. This poses the question that whether there is a system-size threshold for CAC-CIM-CDP in the very-large-scale regime. Considering the MRI-based simulations require 4096 and 16384 DOPO pulses to operate (compared to 16 DOPOs in theoretical simulations in Inui2022 ), the system size of CAC’s applicability is largely improved. Yet the system-size-wise dependency is yet to be explored.

3.2 Advantages of CAC-CIM architecture

With the use of CDP, the problem which involves quadratic optimisation has been solved in this hybrid system. As shown in the schematic illustration of the CAC-CIM-CDP in Fig. 1, proposing approach performs an alternating minimisation between the CIM and CDP. It is clear considering the results stated in Section 2.5 that CAC-CIM-CDP has outperformed OL-CIM-CDP and the generally used approach LASSO which is an $l_{1}$-regularised method for solving compressed sensing problems. It is interesting to see that advancements in CIM architecture can offer better results in real-world problem instances.

3.3 CAC-CIM-CDP (Wigner) Vs. CAC-CIM-CDP (Positive-$P$)

Even though this paper introduces two variants (Wigner and Positive-$P$) of CAC-CIM-CDP, the performances have been almost identical between the models. However, we encountered a deviation when the problem instances become harder i.e. sparseness/compression ratio becomes higher when $w_{noise}=0$. The results are presented in Supplementary note 2. As the models approach a threshold point for optimal reconstruction (a critical sparseness/compression ratio), beyond that the producing RMSE values are somewhat different between the models. Performance-wise it is hard to state that one model is better than the other. Because the significance of Wigner and Positive-$P$ lies in the density matrix approximation and how it behaves with a large quantum noise presence. We discuss this in Supplementary note 2.

3.4 Future improvements to the CAC-CIM-CDP

3.4.1 Simultaneous minimisation

One of the major bottlenecks the proposed model (CAC-CIM-CDP) has is the alternating minimisation process between the CIM and CDP. This is a time-consuming operation. As a future direction to this model, we plan to improvise the CIM system to accommodate quadratic optimisation problems and perform simultaneous minimisation using only the CIM to solve compressed sensing problems. We believe that the use of ”CIM-only” will have a positive effect on accuracy as well.

3.4.2 CAC-CIM-CDP with large quantum noise

While this manuscript solely focuses on combining CAC with CIM for solving CS problems more accurately, the considered quantum noise present in the CIM is very low ($g^{2}=10^{-7}$). This opens up a problem of whether CAC-CIM-CDP can keep up the performance with a large quantum noise presence. For small-scale frustrated Ising hamiltonians, this has been previously explored in Inui2022 ($N=16$) where it has shown a decrease in success probability for larger $g^{2}$ terms. This result is consistent with CAC-CIM-CDP as well as shown in Supplementary note Fig. LABEL:mri64saturationgraphs for MRI simulations. Recent advances in CIM research have led to the introduction of a method known as Negative Parametric Gain (NPG), which accommodates higher quantum noise and at the same time as maintaining a higher probability of success npg . This method considers a negative starting pump rate with large injection field feedback. NPG has shown promising results in the theoretical simulations npg . We are planning to improve the endurance of the CAC-CIM-CDP with NPG for a larger quantum noise presence.

3.4.3 CAC-CIM-CDP with the mean-field CIM model

As it is obvious from the perspective of numerical simulations, CAC-CIM-CDP SDEs are computationally costly to simulate. Even though the shown results are acquired using a GPU implementation of the SDEs, as a digital simulator, field-programmable gate arrays (FPGAs) are more suitable (less energy cost, faster processing etc). As a future direction, we plan on implementing the mean-field CIM SDEs sam ; AmpLeleu with CAC on an FPGA to perform compressed sensing simulations. Due to the fact that CAC-CIM-CDP has relatively low noise present in the system, we believe that the mean-field SDEs will have approximately the same or better results but with faster simulation times. This is mainly due to the simplicity and the negligence of the noise terms in the mean-field CIM SDEs.

4 Methods

4.1 Stochastic Differential equation in OL-CIM-CDP and CAC-CIM-CDP

4.1.1 Wigner-type

The CIM model based on the Wigner formulation was introduced in nonoiseCIM ; Inui . The $c$-number Heisenberg Langevin equation nonoiseCIM was used to overcome the higher computational cost of simulating the direct density matrix formulation of CIM and it has been found to be equivalent to the truncated Wigner SDEs. The density operator master equation expanded by the Wigner function results in the Kramers-Moyal series including third-order terms. In order to derive the Langevin equation, we neglect third-order terms Inui2022 . Then, we can formulate the following Wigner SDEs used for OL-CIM-CDP.

[TABLE]

[TABLE]

Here, in-phase and quadrature-phase normalised amplitudes are represented as $c_{r}$ and $s_{r}$ respectively. $p$ is the normalised pump rate. If $p$ is above the oscillation threshold $(p>1)$, each of the OPO pulses is either in the [math]-phase state or $\pi$-phase state. The last terms of the upper and lower equations express the vacuum fluctuations injected from external reservoirs and the pump fluctuations coupled to the OPO system via gain saturation Aonishi . $W_{r,1}$ and $W_{r,2}$ are independent real Gaussian noise processes satisfying $\langle W_{r,k}(t)\rangle=0$ and $\langle W_{r,k}(t)W_{r^{\prime},l}(t^{\prime})\rangle=\delta_{rr^{\prime}}\delta_{kl}\delta(t-t^{\prime})$. ${g}$ indicates the saturation parameter. $(dc_{r}/dt)_{inj,r}$ is the optical injection field, which only considers the in-phase amplitudes for the calculations. The injection field is defined in Eq. (5) and Eq. (6). $\tilde{K}$ indicates the normalised feedback strength.

Focusing on the behaviour of the OPO pulses only in the in-phase direction, the Wigner-type SDE, which is used for CAC-CIM-CDP, can be stated as,

[TABLE]

[TABLE]

Here $\mu_{r}$ and $V_{r}$ are the mean-amplitudes and the variance of the $r$-th DOPO pulse. $(d\mu_{r}/dt)_{inj,r}$ is the optical injection field defined in Eq. (8)-(11). $W_{R,r}$ is independent real Gaussian noise processes satisfying $\langle W_{R,r}(t)\rangle=0$ and $\langle W_{R,r}(t)W_{R,r^{\prime}}(t^{\prime})\rangle=\delta_{rr^{\prime}}\delta(t-t^{\prime})$. $g$, $p$, $j$ and ${K}$ indicate the saturation parameter, pump rate, the normalised out-coupling rate for optical homodyne measurement and the feedback strength, respectively.

4.1.2 Positive-$P$-type

Positive-$P$ (P-P) representation drumOps is a generalised form of Glauber–Sudarshan $P$ representation. When the density operator master equations are expanded using the P-P distribution function, the resulting Kramers-Moyal series only consists of first and second-order terms. Due to this factor, there is no truncation needed to derive the Langevin equation. Because of this one can argue that P-P SDEs might be a better candidate for density operator approximations. The effectiveness of P-P SDEs has been demonstrated on CIMs with higher quantum noise presence Inui2022 . We can formulate the P-P-type SDEs we used for CAC-CIM-CDP.

[TABLE]

[TABLE]

[TABLE]

Here $\mu_{r}$ corresponds to the mean-amplitude, $m_{r}$ and $n_{r}$ represent variances of quantum fluctuations of the $r$-th DOPO pulse. $(d\mu_{r}/dt)_{inj,r}$ is the optical injection field defined in Eq. (8)-(11). $W_{R,r}$ is independent real Gaussian noise processes satisfying $\langle W_{R,r}(t)\rangle=0$ and $\langle W_{R,r}(t)W_{R,r^{\prime}}(t^{\prime})\rangle=\delta_{rr^{\prime}}\delta(t-t^{\prime})$. $g$,$p$, $j$ and ${K}$ are the same as those for the Wigner model.

4.2 Optimisation in CDP

The CDP performs the optimisation of the Hamiltonian (Eq. 4) with respect to $R_{r}$ for a support vector $\sigma$ given by CIM. $\sigma$ is obtained by binarising the measured-amplitude ($\tilde{\mu}_{r}$) defined in Eq. (10) (CAC-CIM-CDP) or in-phase amplitude $c_{r}$ (OL-CIM-CDP) with the Heaviside function stated as,

[TABLE]

The CDP solve the following system of equations, which is satisfied the stationary point that minimises $\mathbb{H}$ with respect to $r$.

[TABLE]

[TABLE]

Here, $\mathbb{H}_{r}$ in Eq. (22) is the local field of the CDP, which is the same as Eq. (4) and (11). For the simulations, we used the Jacobi method or Conjugate Gradient Descent (CGD) method as the CDP optimiser. During the optimisation in the CDP, all $\sigma_{r}$ are fixed.

4.3 Schedule of pump rate, threshold and target amplitude for optimisation in CIM

A rough parameter search was used to determine the schedules for each of the following parameters in the experiments. The pump rate $p$ for both Wigner and P-P type CAC-CIM-CDPs was scheduled depending on the time $t$ as follows.

[TABLE]

Here, $p_{thr}=1$ for all simulations of both Wigner and P-P type CAC-CIM-CDPs. For artificial random data and MRI data simulations, $d$ was set at 0.6 and 0.4 respectively.

In accordance with Aonishi , the pump rate $p$ for OL-CIM-CDP was scheduled depending on the time $t$ as follows.

[TABLE]

The pump rate becomes equal to $1.5$ when $t=5$. We used this pump rate schedule for all simulations of OL-CIM-CDP. In both CAC-CIM-CDP and OL-CIM-CDP, the threshold $\eta$ was scheduled depending on the alternating iteration time $i$ as follows.

[TABLE]

Here $velo=51$ for all simulations of both CAC-CIM-CDP and OL-CIM-CDP in artificial random data. For the MRI data, $velo=31$ and $velo=11$ were used in OL-CIM-CDP and CAC-CIM-CDP respectively. For synthesised random data (Figs. 3, 4, and Supplementary note Fig. LABEL:GACSsupport), the threshold $\eta$ was linearly lowered from $\eta_{init}$ to $\eta_{end}$ as the alternating minimisation proceeds. $\eta_{init}$ and $\eta_{end}$ are adjusted to maximise the performance of those models. On the other hand, for MRI data (Figs. 5, 6, 7, Supplementary note Fig. LABEL:mri64saturationgraphs and Supplementary note Fig. LABEL:mri128maingraphs), the threshold $\eta$ was constant by setting as $\eta_{init}=\eta_{end}$. The values of $\eta_{init}$ and $\eta_{end}$ used for each simulation are shown in the figure captions.

In both Wigner and P-P type CAC-CIM-CDPs, the target amplitude of CAC, $\tau$, was constant with respect to the time $t$. For the simulations in Fig. 4a and Fig. 4b, $\tau=0.21$ was used while in Fig. 4c and Fig. 4d $\tau$ was set to 0.15. For other simulations, $\tau$ was $1$.

4.4 Observation model for Compressed Sensing

The observation model that is the premise of Eq. (3) and Eq. (4) is defined as follows.

[TABLE]

Here, $A\in\mathbb{R}^{N\times M}$ is the observation matrix, $y\in\mathbb{R}^{M}$ implies the observation signal, $x\in\mathbb{R}^{N}$ and $\xi\in(0,1)^{N}$ are the true source signal and true support, respectively. $w_{noise}\in\mathbb{R}^{M}$ indicates the observation noise satisfying $\langle w_{noise}^{k}\rangle=0$ and $\langle w_{noise}^{k}w_{noise}^{k^{\prime}}\rangle=\nu^{2}\delta_{kk^{\prime}}$. $\nu^{2}$ is the variance of the observation noise.

4.5 Artificial random data

To verify the performance of the proposed models statistically and moreover compare those results with ground states predicted with statistical mechanics Aonishi , we used many samples of artificial random data $y\in\mathbb{R}^{M}$ synthesised from the observation model Eq. (26) in which the values of all entries were randomly determined as follows. Each entry of the observation matrix $A\in\mathbb{R}^{M\times N}$ is randomly generated from an independent and identical normal distribution with the variance of $1/M$, which satisfies $\langle A_{r}^{k}\rangle=0$ and $\langle A_{r}^{k}A_{r^{\prime}}^{t^{\prime}}\rangle=1/M\delta_{rr^{\prime}}\delta(kk^{\prime})$.

Each entry of the true source signal $x\in\mathbb{R}^{N}$ is randomly generated from an independent and identical normal distribution with the variance of $1$, which satisfies $\langle x_{r}\rangle=0$ and $\langle x_{r}x_{r^{\prime}}\rangle=\delta_{rr^{\prime}}$. $a\times N$ elements of $\xi\in(0,1)^{N}$ are randomly selected and assigned $1$ while others are assigned [math]. $a$ is the sparseness defined in the Introduction.

4.6 Simulations with MRI data

To evaluate the performance of the proposed models on realistic data, we used MRI data provided from the fastMRI datasets fastmri . The initial brain MRI used here was a $320\times 320$ image. To reduce the problem size, we resized the image to $64\times 64$ and $128\times 128$ images with the BILINEAR interpolation method. We applied the Haar-wavelet transform (HWT) to the two different-sized images and in Fig. 6 and Fig. 7 we set 78.8% and 82.2% of the HWT coefficients to zero to create two different-sized sparse images ($64\times 64$ and $128\times 128$ pixels) spanned by Haar basis functions with a sparseness of 0.212 and 0.178, respectively. Then, we applied the discrete Fourier transform (DFT) to the two different-sized sparse images to obtain $64\times 64$ and $128\times 128$ $k$-space data, respectively. Finally, we undersampled 1638 and 4915 points from the $64\times 64$ and $128\times 128$ $k$-space data at random red points (Fig. 6b and Fig. 7b) to create two observation signals with a compression rate of 0.4 and 0.3 respectively.

In accordance with our previous work, we sought to reconstruct the source signals from the undersampled $k$-space data by solving the following optimisation problem with CAC-CIM-CDP and OL-CIM-CDP.

[TABLE]

Here, $x$ is a source signal, and $y$ is the observation signal constructed through the above steps. $F$ indicates the DFT matrix and $\Psi$ is the HWT matrix. $F$ and $\Psi$ are orthogonal matrices and their transpose matrices correspond to inverse DFT and inverse HWT, respectively. $S$ is an undersampling matrix executing undersampling at random red points shown in Fig. 6b and Fig. 7b. $\Delta_{v}$ and $\Delta_{h}$ are the matrices discretely representing the vertical and horizontal second-order derivative operators, respectively. $\gamma$ and $\lambda$ are the $l_{2}$ and $l_{0}$ regularisation parameters.

To implement the optimisation problem in Eq. (27) on CIM, we estimate the HWT coefficients instead of the pixel values of the image. Applying the HWT $r=\Psi x$ to Eq. (27), the mutual interaction matrix $J$ and the Zeeman term vector $h^{z}$ for CIM are given as

[TABLE]

[TABLE]

Here, the observation matrix is given as $A=SF\Psi^{T}$. The second and third terms in $\tilde{J}$ are from the $l_{2}$ regularisation terms. After the alternating minimisation, the output of the CDP, $r$, is transformed to the image, $x$, with the inverse HWT $x=\Psi^{T}r$. $\gamma$ is set to 0.0001. $\tilde{K}$ for OL-CIM-CDP was set to 0.25 while ${K}$ for CAC-CIM-CDP was 0.01. Here we use LASSO’s solution as the initial condition for the CIM simulation.

\bmhead

Authors’ contributions

M.D.S.H.G., and T.A., modelled the system, performed the numerical simulations for the proposed models and wrote the manuscript. M.D.S.H.G., T.A., and S.K., worked on the evaluation of the models. K.M., and M.O., provided feedback on numerical simulations. S.K., Y.I., and Y.Y., helped with the physics of CIM and provided feedback.

\bmhead

Availability of data and materials

The data generated and/or analysed during this study are not publicly available for legal/ethical reasons. But M.D.S.H.G. can provide the raw data if formally requested.

\bmhead

Competing Interests

The authors declare no competing interests.

\bmhead

Funding

This work is supported by the Japan Science and Technology Agency through its ImPACT program, NTT Research Inc. And Authors acknowledges the support of the NSF CIM Expedition award (CCF-1918549).