Quantum Reservoir Computing

1. 1
   Step 1:

   $$x_k=F_{s_k}\circ F_{s_{k-1}}\circ\cdots\circ F_{s_{k-K+1}}\big(x_{k-K}\big),\qquad F_s(x)=F(x,s)$$

   Unroll the last $K$ updates. The state at time $k$ depends on the past only through the state $x_{k-K}$ reached $K$ steps earlier.
2. 2
   Step 2:

   $$x'_k=F_{s_k}\circ\cdots\circ F_{s_{k-K+1}}\big(x'_{k-K}\big)$$

   A second history with possibly different early inputs and a different initial state, but the same last $K$ inputs $s_{k-K+1},\dots,s_k$, uses the same $K$ maps applied to another state $x'_{k-K}$.
3. 3
   Step 3:

   $$D(x_k,x'_k)\le\sup_{a,b\in\mathcal X}D\big(F_{s_k}\circ\cdots\circ F_{s_{k-K+1}}(a),\,F_{s_k}\circ\cdots\circ F_{s_{k-K+1}}(b)\big)\le\delta_K$$

   Both are $K$-step trajectories of the same inputs from $a=x_{k-K}$ and $b=x'_{k-K}$; the uniform echo-state property bounds their distance by $\delta_K$ for every pair of initial states.
4. 4
   Step 4:

   $$\delta_K\le\varepsilon\ \text{ for }K\ge K(\varepsilon)\quad\Rightarrow\quad D(x_k,x'_k)\le\varepsilon$$

   Since $\delta_K\to0$, choose the horizon $K(\varepsilon)$ with $\delta_K\le\varepsilon$. This is the fading-memory property, and the feature map turns it into a statement about the features if $\phi$ is continuous.
5. 5
   Step 5:

   $$\big\lvert\langle Z_q\rangle_{\rho}-\langle Z_q\rangle_{\sigma}\big\rvert=\big\lvert\mathrm{Tr}\big[Z_q(\rho-\sigma)\big]\big\rvert\le\lVert Z_q\rVert_\infty\,\lVert\rho-\sigma\rVert_1=2\,T(\rho,\sigma)$$

   For the quantum features: Hölder's inequality with $\lVert Z_q\rVert_\infty=1$ and the trace distance $T=\tfrac12\lVert\rho-\sigma\rVert_1$. Each feature is Lipschitz in the state, so states within $\varepsilon$ give features within $2\varepsilon$.

1. 1
   Step 1:

   $$\hat y_k=\phi_k^\top w,\qquad\phi_k=\phi(x_k)\ \text{ independent of }w$$

   The states, and therefore the feature vectors $\phi_k$, are computed once from the fixed dynamics; the prediction is linear in the trainable parameters.
2. 2
   Step 2:

   $$L(w)=\sum_k\big(\phi_k^\top w-y_k\big)^2+\lambda\lVert w\rVert^2=\lVert Xw-y\rVert^2+\lambda\lVert w\rVert^2,\qquad X_{kj}=(\phi_k)_j$$

   The squared-error loss with an $\ell_2$ penalty is a quadratic polynomial in $w$; stacking the feature vectors as the rows of a matrix $X$ gives the compact form.
3. 3
   Step 3:

   $$\nabla^2L=2\big(X^\top X+\lambda I\big)\succ0$$

   The Hessian is constant and positive definite for $\lambda>0$ (ridge section), so $L$ is strictly convex: no local minima, no learning rate, and a closed-form minimiser.
4. 4
   Step 4:

   $$\frac{\partial x_k}{\partial\theta}=\sum_{m\le k}\Big(\prod_{j=m+1}^{k}\frac{\partial F}{\partial x}\Big|_{x_{j-1},s_j}\Big)\frac{\partial F}{\partial\theta}\Big|_{x_{m-1},s_m}$$

   Training a parameter $\theta$ of the dynamics instead would require this backpropagation-through-time sum, with products of up to $k$ Jacobians.
5. 5
   Step 5:

   $$\Big\lVert\prod_{j=m+1}^{k}\frac{\partial F}{\partial x}\Big\rVert\le\prod_{j=m+1}^{k}\eta_j\ \xrightarrow{\ k-m\to\infty\ }\ 0$$

   For the quantum reservoir each step is a linear map on density matrices whose trace-norm gain $\eta_j$ on differences of states is at most $1$ (next section). The contraction that produces fading memory also makes these gradient products vanish, and on hardware the state $x_k$ is not even accessible.

1. 1
   Step 1:

   $$H=\sum_{0\le i<j\le3}J_{ij}\,X_iX_j+h\sum_{q=0}^{3}Z_q$$

   The blueprint model: an Ising interaction in the $x$ direction between every pair of qubits and a field $h$ in the $z$ direction, which does not commute with the couplings and therefore generates nontrivial dynamics.
2. 2
   Step 2:

   $$(J_{01},J_{02},J_{03},J_{12},J_{13},J_{23})=(0.0695,\ 0.6305,\ -0.0920,\ -0.0749,\ -0.3440,\ 0.8774),\qquad h=0.5$$

   The fixed couplings of the lab, drawn once at random and then frozen. Couplings of different sizes and signs avoid symmetries between qubits, which would make some features linearly dependent.
3. 3
   Step 3:

   $$H=0.0695\,XXII+0.6305\,XIXI-0.0920\,XIIX-0.0749\,IXXI-0.3440\,IXIX+0.8774\,IIXX+0.5\,(ZIII+IZII+IIZI+IIIZ)$$

   The same operator as a sum of ten Pauli words, which is how the engine stores it before building the $16\times16$ matrix.
4. 4
   Step 4:

   $$H^\top=H=\overline H,\qquad\operatorname{spec}H\subset[-2.6358,\ 2.5890]$$

   $X$ and $Z$ are real symmetric, hence so are their tensor products and $H$. Its sixteen eigenvalues lie between $-2.6358$ and $2.5890$ and are non-degenerate. The Bohr frequencies $\omega_{ab}=E_a-E_b$ that appear in $\langle Z_q\rangle(t)$ range from $0.0548$ to $5.2248$, i.e. periods $2\pi/\omega_{ab}$ from $1.203$ to $114.7$: during one input period $\tau=4$ the fastest modes complete more than three oscillations while the slowest barely move, so the features mix fast and slow time scales.
5. 5
   Step 5:

   $$X_iX_j\,P=P\,X_iX_j,\quad Z_q\,P=P\,Z_q,\qquad P=Z_0Z_1Z_2Z_3\ \Rightarrow\ [H,P]=0$$

   $X$ anticommutes with $Z$ on the same qubit, so $X_iX_j$ anticommutes with $Z_i$ and with $Z_j$, and two sign changes cancel. The unitary dynamics alone conserves the parity $P$; the injection, which resets qubit $0$, breaks this conservation law and lets the input steer the state.

1. 1
   Step 1:

   $$\mathcal R(\rho)=\sum_{b}\big(\lvert0\rangle\langle b\rvert_0\otimes I\big)\rho\big(\lvert b\rangle\langle0\rvert_0\otimes I\big)=\lvert0\rangle\langle0\rvert_0\otimes\mathrm{Tr}_0\,\rho$$

   The non-selective reset of qubit $0$ has the Kraus operators $\lvert0\rangle\langle0\rvert$ and $\lvert0\rangle\langle1\rvert$; it is the case $s=0$ of the theorem, since $\lvert\psi_0\rangle=\lvert0\rangle$.
2. 2
   Step 2:

   $$R_Y(\theta)=\begin{pmatrix}\cos\frac\theta2&-\sin\frac\theta2\\\sin\frac\theta2&\cos\frac\theta2\end{pmatrix},\qquad R_Y(\theta)\lvert0\rangle=\cos\tfrac\theta2\lvert0\rangle+\sin\tfrac\theta2\lvert1\rangle$$

   $R_Y(\theta)=e^{-i\theta Y/2}=\cos\frac\theta2I-i\sin\frac\theta2Y$ with $-iY=\begin{pmatrix}0&-1\\1&0\end{pmatrix}$; its first column is the image of $\lvert0\rangle$.
3. 3
   Step 3:

   $$\theta_s=2\arcsin\sqrt s:\qquad\sin\tfrac{\theta_s}2=\sqrt s,\qquad\cos\tfrac{\theta_s}2=\sqrt{1-\sin^2\tfrac{\theta_s}2}=\sqrt{1-s}$$

   $\arcsin\sqrt s\in[0,\frac\pi2]$ for $s\in[0,1]$, where the cosine is non-negative, so the positive square root is the right one. Hence $R_Y(\theta_s)\lvert0\rangle=\lvert\psi_s\rangle$.
4. 4
   Step 4:

   $$R_Y(\theta_s)=\begin{pmatrix}\sqrt{1-s}&-\sqrt s\\\sqrt s&\sqrt{1-s}\end{pmatrix}$$

   The rotation written with the input itself, exactly the matrix the engine applies after its reset.
5. 5
   Step 5:

   $$\big(R_Y(\theta_s)\otimes I\big)\big(\lvert0\rangle\langle0\rvert\otimes\sigma\big)\big(R_Y(\theta_s)\otimes I\big)^\dagger=\lvert\psi_s\rangle\langle\psi_s\rvert\otimes\sigma\ \ \Rightarrow\ \ \mathcal T_s=\mathrm{Ad}_{R_Y(\theta_s)}\circ\mathcal R$$

   Apply the rotation to the reset state with $\sigma=\mathrm{Tr}_0\rho$. The injection is a reset followed by an ordinary data-dependent rotation, as in angle encoding, except that the reset discards whatever qubit $0$ carried.
6. 6
   Step 6:

   $$\langle Z_0\rangle=\mathrm{Tr}\big[(Z\otimes I)\,(\lvert\psi_s\rangle\langle\psi_s\rvert\otimes\sigma)\big]=\langle\psi_s\vert Z\vert\psi_s\rangle\,\mathrm{Tr}\,\sigma=(1-s)-s=1-2s$$

   Right after the injection, the magnetisation of qubit $0$ encodes the input linearly; also $\langle X_0\rangle=2\sqrt{s(1-s)}$ and $\langle Y_0\rangle=0$, so the Bloch vector lies in the $xz$ plane at polar angle $\theta_s$. For the first input $s_0=0.6011$ of the lab, $\langle Z_0\rangle=-0.2022$.

Quick check

1. 1
   Step 1:

   $$U=e^{-iH\tau/V},\qquad U^V=e^{-iH\tau}$$

   All factors are exponentials of the same operator, so they commute and their exponents add. The $V$ sub-steps reproduce the evolution over the full time $\tau$ between inputs; here $\tau/V=0.8$.
2. 2
   Step 2:

   $$\rho_k^{(0)}=\mathcal T_{s_k}(\rho_{k-1}),\qquad\rho_k^{(v)}=U\rho_k^{(v-1)}U^\dagger=U^v\rho_k^{(0)}U^{-v}\quad(v=1,\dots,V),\qquad\rho_k=\rho_k^{(V)}$$

   Inject the input, then evolve in $V$ sub-steps. The state after the last sub-step is the reservoir state that the next injection acts on.
3. 3
   Step 3:

   $$Z_q\lvert x\rangle=(-1)^{x_q}\lvert x\rangle,\qquad x_q=\big\lfloor x/2^{\,n-1-q}\big\rfloor\bmod2$$

   $Z_q$ is diagonal in the computational basis; with qubit $0$ as the most significant bit, the bit of qubit $q$ is bit $n-1-q$ of the integer $x$.
4. 4
   Step 4:

   $$\langle Z_q\rangle_\rho=\mathrm{Tr}(Z_q\rho)=\sum_x\langle x\vert Z_q\rho\vert x\rangle=\sum_x(-1)^{x_q}\rho_{xx}$$

   Evaluate the trace in the computational basis and use $\langle x\rvert Z_q=(-1)^{x_q}\langle x\rvert$. Only the diagonal, i.e. the measurement probabilities $\rho_{xx}$, enters.
5. 5
   Step 5:

   $$\phi_k=\Big(\langle Z_0\rangle_{\rho_k^{(1)}},\dots,\langle Z_3\rangle_{\rho_k^{(1)}},\ \dots,\ \langle Z_0\rangle_{\rho_k^{(V)}},\dots,\langle Z_3\rangle_{\rho_k^{(V)}},\ 1\Big)\in\mathbb R^{nV+1}$$

   The feature row of time step $k$, ordered sub-step by sub-step, followed by a constant $1$ that lets the readout add an offset (the bias column).
6. 6
   Step 6:

   $$nV+1=4\cdot5+1=21,\qquad X\in\mathbb R^{600\times21}\ \text{(training rows }k=200,\dots,799)$$

   The lab readout therefore has $21$ weights. The first $200$ rows are discarded as washout and the last $400$ rows are the test window.
7. 7
   Step 7:

   $$2V\cdot16^3=40\,960\ \text{complex multiply-adds per input}$$

   The cost of the exact simulation: each sub-step computes $U\rho U^\dagger$, two products of $16\times16$ matrices. On hardware each expectation value needs many shots and measurement disturbs the state, so each feature needs an ensemble of identical systems (the NMR setting of Fujii and Nakajima) or a re-run of the input history.

Quick check

1. 1
   Step 1:

   $$\sum_k\big(\hat y_k-y_k\big)=-\lambda\,W_{F-1}$$

   Row $F-1$ (the bias column, all ones) of the normal equations reads $\sum_kX_{k,F-1}\big((XW)_k-y_k\big)+\lambda W_{F-1}=0$. Because the bias sits inside $X$, the penalty shrinks it as well, as in the engine. Textbook ridge often leaves the intercept unpenalised, which would make the mean residual exactly zero.
2. 2
   Step 2:

   $$\frac1T\sum_k(\hat y_k-y_k)=-\frac{10^{-6}\cdot0.4913}{600}=-8.19\times10^{-10}$$

   For the lab's delay-$1$ readout ($\lambda=10^{-6}$, $T=600$, bias weight $0.4913$) the effect is negligible: the mean residual on the training rows is below $10^{-9}$.
3. 3
   Step 3:

   $$A=X^\top X+\lambda I=LL^\top,\qquad A_{ij}=\sum_{m\le j}L_{im}L_{jm}\quad(i\ge j)$$

   Every symmetric positive definite matrix has a unique Cholesky factorisation with $L$ lower triangular and $L_{jj}>0$. Writing out the product shows which entries of $L$ enter each $A_{ij}$.
4. 4
   Step 4:

   $$L_{jj}=\sqrt{A_{jj}-\sum_{m<j}L_{jm}^2},\qquad L_{ij}=\frac1{L_{jj}}\Big(A_{ij}-\sum_{m<j}L_{im}L_{jm}\Big)\quad(i>j)$$

   Solve the previous equation column by column for the diagonal entry and then for the entries below it. The radicand is positive: it equals $\det A_{[0..j]}/\det A_{[0..j-1]}$, a ratio of leading principal minors, which are positive for a positive definite matrix.
5. 5
   Step 5:

   $$Lz=X^\top y:\ z_i=\frac{(X^\top y)_i-\sum_{m<i}L_{im}z_m}{L_{ii}};\qquad L^\top W=z:\ W_i=\frac{z_i-\sum_{m>i}L_{mi}W_m}{L_{ii}}$$

   $AW=L(L^\top W)$, so the system splits into a forward substitution from the top and a back substitution from the bottom.
6. 6
   Step 6:

   $$TF^2+\tfrac13F^3+2F^2=600\cdot441+3087+882\approx2.7\times10^5\ \text{operations}$$

   Forming $X^\top X$ dominates. The factorisation costs $F^3/3$ multiply-adds and the two triangular solves $F^2$ each; the lab's system is tiny.
7. 7
   Step 7:

   $$\operatorname{spec}\big(X^\top X\big)\subset[0.002504,\ 614.7]\ \Rightarrow\ \kappa\approx2.5\times10^5$$

   The measured eigenvalues of the lab's $21\times21$ Gram matrix. The condition number is large but harmless in double precision; $\lambda=10^{-6}$ shifts the smallest eigenvalue by less than $0.1\%$, so it guards against singularity without biasing the fit.

1. 1
   Step 1:

   $$\hat y_k\equiv c:\qquad\sum_k(c-y_k)^2=\sum_k(y_k-\bar y)^2+T(c-\bar y)^2$$

   Insert $c-y_k=(c-\bar y)-(y_k-\bar y)$ and expand; the cross term $2(c-\bar y)\sum_k(y_k-\bar y)$ vanishes by the definition of $\bar y$.
2. 2
   Step 2:

   $$\mathrm{NMSE}_{c}=1+\frac{(c-\bar y)^2}{\mathrm{Var}(y)},\qquad\mathrm{NMSE}_{\bar y}=1$$

   Divide by $\sum_k(y_k-\bar y)^2=T\,\mathrm{Var}(y)$. Predicting the mean scores exactly $1$; a useful readout must stay well below $1$.
3. 3
   Step 3:

   $$\hat y\equiv0:\ \mathrm{NMSE}_0=1+\frac{\bar y^2}{\mathrm{Var}(y)}\approx1+\frac{1/4}{1/12}=4\quad(s\sim\mathcal U[0,1))$$

   The all-zero readout. For inputs uniform on $[0,1)$ the mean is $\frac12$ and the variance $\frac1{12}$. On the lab's test window $\bar y=0.5131$ and $\mathrm{Var}(y)=0.08075$, which gives $4.260$.
4. 4
   Step 4:

   $$\hat y_k=s_{k-d'},\ y_k=s_{k-d},\ d'\ne d:\qquad\frac1T\sum_k(s_{k-d'}-s_{k-d})^2\approx2\,\mathrm{Var}(s)\ \Rightarrow\ \mathrm{NMSE}\approx2$$

   Recalling the wrong delay perfectly is worse than predicting the mean: for independent inputs of equal variance, $\mathbb E(a-b)^2=\mathrm{Var}\,a+\mathrm{Var}\,b+(\mathbb Ea-\mathbb Eb)^2=2\,\mathrm{Var}(s)$.
5. 5
   Step 5:

   $$\lambda=0:\ X^\top(y-\hat y)=0\ \Rightarrow\ \sum_k(y_k-\hat y_k)=0,\quad\sum_k(\hat y_k-\bar y)(y_k-\hat y_k)=0$$

   On the training rows of an ordinary least-squares fit, the residual is orthogonal to every column of $X$. The bias column gives a zero mean residual, and since $\hat y-\bar y\mathbf 1=X(W-\bar ye_{F-1})$ lies in the column space, it is orthogonal to the residual.
6. 6
   Step 6:

   $$\sum_k(y_k-\bar y)^2=\sum_k(\hat y_k-\bar y)^2+\sum_k(y_k-\hat y_k)^2\ \Rightarrow\ r^2(\hat y,y)=\frac{\sum_k(\hat y_k-\bar y)^2}{\sum_k(y_k-\bar y)^2}=1-\mathrm{NMSE}$$

   Pythagoras for the orthogonal decomposition $y-\bar y=(\hat y-\bar y)+(y-\hat y)$. Moreover $\sum_k(\hat y_k-\bar y)(y_k-\bar y)=\sum_k(\hat y_k-\bar y)^2$, so the squared Pearson correlation equals the explained fraction. On the training rows $r^2=1-\mathrm{NMSE}$ exactly; on test rows only approximately.