Have you tryed a simple newton algorithm?

Let $(\alpha_{k})$ be defined as: $\alpha_k=1/k^2$

Initialisation :

$x^0=[0,...,0]$

Compute $H=(A^* A)^{-1}$

Loop for $k$ in $1:m$

$x^k=x^{k-1}-\alpha_k H A^*(Ax^{k-1}-b)$

$x^{k}=g*x^{k}/\|x_2^k\|$

end for loop

Obviously, there are more adaptive way of choosing $\alpha_k$... but maybe you don't need such sofistication to solve a norm minimization problem. If $A^* A$ has very small eigen values you can use $H_k=(A^* A+\epsilon_k)^{-1}$ instead of $H$ ($\epsilon_k$ decreasing to zero)...

Here is the corresponding R code:

A=t(array(1:1000,c(10,100)))
m=100; b=1:10; g=3; l=5; p=10;
alpha=1:m
alpha=1/alpha^2
x=array(0,c(m,p))
H=t(A)%*%A
svdH=svd(H)
H=svdH$v%*%diag(1/svdH$d)%*%t(svdH$u) 
for (k in 2:m)
{
    x[k,]=x[k-1,]-alpha[k]*H%*%(t(A)%*%(A%*%x[k-1,]-b))
    x[k,]=g*x[k,]/sqrt(sum(x[k,(l+1):p]^2))
    print(sum((A%*%x[k,]-b)^2))
}

Have you tryed a simple newton algorithm (with the constraint added to the algo)?

Let $(\alpha_{k})$ be defined as: $\alpha_k=1/k^2$

Initialisation :

$x^0=[0,...,0]$

Compute $H=(A^* A)^{-1}$

Loop for $k$ in $1:m$

$x^k=x^{k-1}-\alpha_k H A^*(Ax^{k-1}-b)$

$x^{k}=g*x^{k}/\|x_2^k\|$

end for loop

Obviously, there are more adaptive way of choosing $\alpha_k$... but maybe you don't need such sofistication to solve a norm minimization problem. If $A^* A$ has very small eigen values you can use $H_k=(A^* A+\epsilon_k)^{-1}$ instead of $H$ ($\epsilon_k$ decreasing to zero)...

Note that this type of code is relatively general when you want to find the saddle point of a lagrangian and you know how to find the maxima with respect to Lagrange multipliers (in the dual space) (second step of the loop) but you need a gradient descent (or here Newton algo) for finding the minima in the principal space.

Here is the corresponding R code:

A=t(array(1:1000,c(10,100)))
m=100; b=1:10; g=3; l=5; p=10;
alpha=1:m
alpha=1/alpha^2
x=array(0,c(m,p))
H=t(A)%*%A
svdH=svd(H)
H=svdH$v%*%diag(1/svdH$d)%*%t(svdH$u) 
for (k in 2:m)
{
    x[k,]=x[k-1,]-alpha[k]*H%*%(t(A)%*%(A%*%x[k-1,]-b))
    x[k,]=g*x[k,]/sqrt(sum(x[k,(l+1):p]^2))
    print(sum((A%*%x[k,]-b)^2))
}