Optimal threshold decay for L1 homotopy basis pursuit

To solve an underdetermined system of equations where we believe the optimal solution should be the one with minimal L1 norm, we solve the basis pursuit problem: $$\min_{\mathbf{x}}{|\mathbf{x}|_1}$$ s.t. $$\mathbf{Ax}=\mathbf{b}$$.

In practice, this is effectively the same as solving the basis pursuit denoising/LASSO problem $$\min_{\mathbf{x}}{\frac{1}{2}||\mathbf{Ax}-\mathbf{b}||_2^2 + \lambda|\mathbf{x}|_1}$$ with a small $$\lambda$$.

Many methods have been proposed to solve such problems. For very large systems (tens or even hundreds of billions of nonzero components in $$\mathbf{x}$$, solved on a computer cluster), where $$\mathbf{A}$$ and its adjoint are not available as explicit matrices, and where $$\mathbf{b}$$ is real but $$\mathbf{A}$$ and $$\mathbf{x}$$ are complex, I have found an L1 homotopy continuation approach to be reasonably successful, where I solve the LASSO problem using ISTA with a threshold ($$\lambda$$) that decays each iteration, from $$|\mathbf{A}^\dagger\mathbf{b}|_\infty$$ in the first iteration to 0 in the last. I generally get a reasonable solution if I decay exponentially over several hundred iterations, but the minimum number of iterations that is necessary seems to be problem dependent.

Is there any practical theory on what the optimal way to decay the threshold is? Or are there other solver methods that might be more suitable for such large systems?