Implementing minimization in SciPy

Question

I am trying to implement the 'Iterative hessian Sketch' algorithm from https://arxiv.org/abs/1411.0347 page 12. However, I am struggling with step two which needs to minimize the matrix-vector function.

Imports and basic data generating function

import numpy as np
import scipy as sp
from sklearn.datasets import make_regression
from scipy.optimize import minimize
import matplotlib.pyplot as plt
%matplotlib inline

from numpy.linalg import norm

def generate_data(nsamples, nfeatures, variance=1):
    '''Generates a data matrix of size (nsamples, nfeatures)
    which defines a linear relationship on the variables.'''
    X, y = make_regression(n_samples=nsamples, n_features=nfeatures,\
                        n_informative=nfeatures,noise=variance)
    X[:,0] = np.ones(shape=(nsamples)) # add bias terms
    return X, y

To minimize the matrix-vector function, I have tried implementing a function which computes the quanity I would like to minimise:

def f2min(x, data, target, offset):
    A = data
    S = np.eye(A.shape[0])
    #S = gaussian_sketch(nrows=A.shape[0]//2, ncols=A.shape[0] )
    y = target
    xt = np.ravel(offset)


    norm_val = (1/2*S.shape[0])*norm(S@A@(x-xt))**2
    #inner_prod = (y - A@xt).T@A@x

    return norm_val - inner_prod

I would eventually like to replace S with some random matrices which can reduce the dimensionality of the problem, however, first I need to be confident that this optimisation method is working.

def grad_f2min(x, data, target, offset):
A = data
y = target
S = np.eye(A.shape[0])
xt = np.ravel(offset)
S_A = S@A
grad = (1/S.shape[0])*S_A.T@S_A@(x-xt) - A.T@(y-A@xt)

return grad

    x0 = np.zeros((X.shape[0],1))
xt = np.zeros((2,1))
x_new = np.zeros((2,1))
for it in range(1):
    result = minimize(f2min, x0=xt,args=(X,y,x_new),
         method='CG', jac=False )
    print(result)
    x_new = result.x

I don't think that this loop is correct at all because at the very least there should be some local convergence before moving on to the next step. The output is:

    fun: 0.0
     jac: array([ 0.00745058,  0.00774882])
 message: 'Desired error not necessarily achieved due to precision loss.'
    nfev: 416
     nit: 0
    njev: 101
  status: 2
 success: False
       x: array([ 0.,  0.])

Does anyone have an idea if:

(1) Why I'm not achieving convergence at each step

(2) I can implement step 2 in a better way?