indicator function in objective function with $L_2$ norm I am trying to solve an optimization problem. The objective function is as follows:
$\arg \min \lVert\mathbb{A}\mathbf{x} - \mathbf{b}\rVert^2 + \text{other linear least squares terms} + \mathcal{I}(\mathit{x_0<a}) \lVert\mathit{x_0 - a}\rVert^2 +  \mathcal{I}(\mathit {x_n>b}) \lVert\mathit{x_n-b}\rVert^2$
where $\mathcal{I}$is the indicator function that returns $1$ for true condition and $0$ otherwise.
$x_0,x_1,...,x_n$ should be between $a$ and $b$. If $x_0$ or $x_n$ is out of the range, one cost will be added to the objective function.
If the indicator function doesn't appear in the objective function, it's simply one linear least squares optimization problem and is simple to solve. Indicator function is not a continuous function and makes the problem difficult.
I am not an expert on numerical optimization. I search on the internet and it seems that indicator function is used often in deep learning. Any hints, links and materials are appreciated.
 A: Relax (as in perform relaxation of your initial constraints). Consider reformulating your the cost function as having two ReLU components one regarding $(x_0, a)$ and another for $(x_n, b)$:
$$
arg\min \lVert\mathbb{A}\mathbf{x} - \mathbf{b}\rVert^2 + \lambda_a\text{ReLU}(a-x_0) + \lambda_b\text{ReLU}(x_n-b)
$$
where through hyper-parameters $\lambda_a$ and $\lambda_b$ we can control the penalisation for being out of the bound $a$ and $b$ respectively. Pick  $\lambda_a$, $\lambda_b$ to be on the same magnitude as your expected loss values for starters and if you get good results, consider amping them to constraint your optimisation space even more aggressively. There are more formal ways of setting $\lambda$ (e.g. see Cotter et al. (2019) "Two-Player Games for Efficient Non-Convex Constrained Optimization" but they are a bit of an overkill to start of).
In addition to the above and as totally separate take, we can consider using a derivative-free optimization method altogether. We can start with something basic like a a simulated annealing/random search approach. I would suggest initialising these search with the OLS solution $x_{\text{init}}$ that is found by treating the problem as having box-constraints on $[a,b]$ and then feed that  $x_{\text{init}}$ to the simulated annealing (the initial solution can be through  L-BFGS-B or any other "off-the-shelf" constrained optimisation algorithm). The NLOpt project has a variety of algorithms to consider aside simulated annealing.
In general, if you have an optimisation problem that does not have an obvious solution, either approximate it (e.g. with the ReLU components mentioned first), or work around it altogether (e.g. skip the use of gradient information with derivative-free optimization methods). (A third option would be to reformulate it e.g. by using an Augmented Lagrangian approach but that's yet another game.)
A: I see two possible algorithms to solve this problem.
R's Glmnet: Lasso and elastic-net regularized generalized linear models"
Glmnet is a package for R to perform elastic net regularisation. With this package you can also constrain parameters by defining a minimum or maximum for parameters and you can in addition choose different scaling of the penalty term for different coefficients (e.g in your problem only two coefficients are penalised).
Then you can augment you matrix $\mathbb{A}$ with copies of the 0-th and n-th column, and for these you have unpenalized coefficients $x_0^\prime$ and $x_n^\prime$. These coefficients will be constrained. That is: the coefficient $x_0^\prime$ can't exceed $a$ and the coefficient $x_n^\prime$ must be below $b$.
These copies will not create an overdetermined system. The $x_0^\prime$ and $x_n^\prime$ are constrained and the $x_0$ and $x_n$ are penalised.
What will happen with this minimization is that the $x_0$ and $x_n$ remain zero as long the $x_0^\prime$ and $x_n^\prime$ are within the boundaries. So in this way you have the penalty 'kick in' automatically when $x_0 + x_0^\prime > a$ or $x_n + x_n^\prime < b$.
Solving four cases
If you feel uncomfortable about constrained regression algorithms (like glmnet uses but there are probably also others packages) then you can do it for your simple case manually (there is only four options based on whether $x_0 < a$ and $x_n > b$).
Now you do not add an additional copy of a column from the matrix but instead subtract a multiple times ($a$ or $b$) from $\mathbb{b}$. After that compare which of the four is the best result. (The four combinations result from either subtracting the column and penalised, or not subtracting and not penalising)
So you compute either when the indicator is 'on' or 'off'. (and when with the indicator turned on you still get a coefficient within the boundaries then it will be even a better cost with the indicator turned off so this comparison will give you the argument with the lowest cost) With the indicator turned off you need to take care that you only count the case when the argument is within the bounds.
A: This problem is quickly and easily solvable if you split up the regions into four pieces:

*

*$x_0\geq a \cap x_n\leq b$: Solve $\arg\min ||\mathbb{A}x-b||^2 + \text{other linear least squares terms}$

*$x_0<a \cap x_n\leq b$: Solve $\arg\min ||\mathbb{A}x-b||^2 + \text{other linear least squares terms} + ||x_0-a||^2$

*$x_0\geq a \cap x_n>b$: Solve $\arg\min ||\mathbb{A}x-b||^2 + \text{other linear least squares terms} + ||x_n-b||^2$

*$x_0<a \cap x_n>b$: Solve $\arg\min ||\mathbb{A}x-b||^2 + \text{other linear least squares terms} + ||x_0-a||^2 + ||x_n-b||^2$
(It is unclear how $x_2,\ldots,x_{n-1}$ are constrained, apart from lying in $[a,b]$.)
The indicator functions are not continuous, but they are multiplying quadratics so the results are still continuous functions. Continuity is not your major concern.
The bigger issue is the second derivatives of those terms are only piecewise-continuous. This will lead methods which approximate the inverse Hessian (like BFGS) to take long to converge (or possibly diverge with a poor starting point). True: many methods could be used. You can do relaxations or you could approximate the indicators with quartics or exponentials or logits. The problem is that turning hard constraints into soft constraints can lead to ill-posed problems. (Otherwise, integer programming would not be NP-hard.)
For that reason, I would recommend splitting the problem up. This is one of the situations where using your understanding to setup the problem can be better than just pressing a button.
