Train the ridge readout of a quantum reservoir

Complete the TypeScript function

trainReadout(features: number[][], targets: number[], lambda: number): number[]

The reservoir (fixed). Four qubits (qubit $0$ is the most significant bit) evolve under $H=\sum_{i<j}J_{ij}X_iX_j+h\sum_qZ_q$ with $J_{01}=0.0695$, $J_{02}=0.6305$, $J_{03}=-0.0920$, $J_{12}=-0.0749$, $J_{13}=-0.3440$, $J_{23}=0.8774$ and $h=0.5$, starting in $\lvert0000\rangle$. At time step $k$ the input $s_k\in[0,1)$ replaces qubit $0$: $\rho\mapsto\lvert\psi_{s_k}\rangle\langle\psi_{s_k}\rvert\otimes\mathrm{Tr}_0\,\rho$ with $\lvert\psi_s\rangle=\sqrt{1-s}\,\lvert0\rangle+\sqrt s\,\lvert1\rangle$. The state then evolves for $\tau=4$ in $V=5$ sub-steps $U=e^{-iH\tau/V}$, and after each sub-step the four values $\langle Z_0\rangle,\dots,\langle Z_3\rangle$ are recorded. The inputs are $s_k=$ mulberry32(42)() for $k=0,\dots,1199$.

The arguments. features is the matrix $X$ with one row per time step $k=200,\dots,799$ (the first $200$ steps are the washout, which removes the dependence on the initial state). Each row has $21$ entries: the $4\times5=20$ expectation values, ordered sub-step by sub-step, followed by a bias entry that is always $1$. Do not add another bias column. targets is the vector $y$ with $y_k=s_{k-d}$, the input $d$ steps in the past, for the same rows in the same order. lambda is the ridge parameter $\lambda=10^{-6}$.

What to compute. The weights $w=\arg\min_{w\in\mathbb R^{21}}\ \lVert Xw-y\rVert^2+\lambda\lVert w\rVert^2 .$ Setting the gradient $2X^\top(Xw-y)+2\lambda w$ to zero gives the normal equations $\big(X^\top X+\lambda I\big)\,w=X^\top y,\qquad (X^\top X)_{ij}=\sum_kX_{ki}X_{kj},\quad (X^\top y)_i=\sum_kX_{ki}\,y_k .$ The $21\times21$ matrix $A=X^\top X+\lambda I$ is symmetric positive definite, so the solution is unique. Solve the system with a Cholesky factorisation $A=LL^\top$ followed by forward substitution $Lz=X^\top y$ and back substitution $L^\top w=z$, or with Gaussian elimination. The regulariser applies to every weight, the bias included. Return $w$ as an array of $21$ finite numbers.

Grading. Your function runs twice, with the delay-$1$ targets (code test weights-d1) and with the delay-$2$ targets (weights-d2); you are not told which delay a call belongs to. The check shape verifies that each call returns $21$ finite numbers. The checks delay-1 and delay-2 recompute the reservoir on the server, predict $\hat y_k=\sum_jX_{kj}w_j$ on the test steps $k=800,\dots,1199$, which your function never sees, and compute $\mathrm{NMSE}=\sum_k(\hat y_k-y_k)^2/\sum_k(y_k-\bar y)^2$. They require $\mathrm{NMSE}\le0.03$ for delay $1$ and $\le0.08$ for delay $2$. The exact ridge solution reaches $0.0105$ and $0.0357$.

About the shortcut. The module 'quantumcore/qml' exports a ridgeRegression function, which the lesson uses. Calling it would pass this lab without doing the exercise. The grader only sees the returned weights, so this rule is not enforced; the point of the lab is to write the normal equations and the solver yourself.