Lagrange-Newton minimization and initial value for the Lagrange multipliers

Short version. How to choose the initial values for the Lagrange multipliers in the Lagrange-Newton equality-constraint minimization method?

Introduction. The problem to solve is

$$\min_{\mathbf x \in \mathbb R^n} f(\mathbf x) \quad \text{subject to} \quad h(\mathbf x) = \mathbf 0,$$

where $f: \mathbb R^n \to \mathbb R$ and $h: \mathbb R^n \to \mathbb R^m$ are twice continuously differentiable functions.

I implemented the local Lagrange-Newton method, which aims at finding a solution $(\mathbf x, \pmb \mu) \in \mathbb{R}^{n+m}$ of the KKT equations

\begin{align} \nabla f(\mathbf x) + \nabla h(\mathbf x)\pmb \mu &= \mathbf 0,\\ h(\mathbf x)&= \mathbf 0, \end{align}

where $\nabla h(\mathbf x) = (\nabla h_1(\mathbf x), \ldots, \nabla h_m(\mathbf x)) \in \mathbb{R}^{n \times m}$ is the Jacobi matrix of $h$ and $\pmb \mu \in \mathbb{R}^m$ is the vector of Lagrange multipliers.

I need initial values $\mathbf x_0$ and $\pmb \mu_0$ in order to use Newton's method to iteratively solve the KKT equations. Based on my application, I have a reasonable value for $\mathbf x_0$.

Question 1 What are good general strategies for selecting $\pmb \mu_0$?

Question 2 How about using $\mathbf x_0$ to compute $\pmb \mu_0$ as the least-squares solution of the first KKT equation: $$\pmb \mu_0 = \text{arg}\min_{\pmb \mu \in \mathbb{R}^m} \|\nabla f(\mathbf x_0) + \nabla h(\mathbf x_0)\pmb \mu\|^2.$$

For some values for $\mathbf x_0$, the least-squares solution $\pmb \mu_0$ makes the method to diverge, whereas setting $\pmb \mu_0$ to a random vector in $[-1, 1]^m$ leads to convergence to the right solution. Any thoughts on that?

Question 3 Is it beneficial to have $\mathbf x_0$ to fulfill the constraints? Or is the distance to the solution $(\mathbf x, \pmb \mu)$ more important?

I think there is an answer to this question in Izmailov AF, Solodov MV: Newton-Type Methods for Optimization and Variational Problems, Springer 2014, section 4.1.1. Newton-Lagrange Method, p. 208.

Here is my attempt at a derivation (in the book, the result is just given). Given the Lagrangian

$$L(\mathbf{x}, \mathbf{\lambda}) = f(\mathbf{x}) + \boldsymbol{\lambda}^T \mathbf{h}(\mathbf{x})$$

we have

$$\frac{\partial L}{\partial\mathbf{x}} = \frac{\partial f}{\partial\mathbf{x}} + \boldsymbol{\lambda}^T \frac{\partial\mathbf{h}}{\partial\mathbf{x}}$$

or

$${\frac{\partial L}{\partial\mathbf{x}}}^T = {\frac{\partial f}{\partial\mathbf{x}}}^T + {\frac{\partial\mathbf{h}}{\partial\mathbf{x}}}^T \boldsymbol{\lambda}.$$

In an extremum this is zero, so we have

$$- {\frac{\partial f}{\partial\mathbf{x}}}^T = {\frac{\partial\mathbf{h}}{\partial\mathbf{x}}}^T \boldsymbol{\lambda}.$$

In the book it is assumed that

$$\frac{\partial\mathbf{h}}{\partial\mathbf{x}}$$

has full column rank, so

$$\frac{\partial\mathbf{h}}{\partial\mathbf{x}} {\frac{\partial\mathbf{h}}{\partial\mathbf{x}}}^T$$

is non-singular.

If we multiply the last equation from the left by

$$\frac{\partial\mathbf{h}}{\partial\mathbf{x}}$$

(which, however, is not an invertible step, so I'm not sure here) we get

$$- \frac{\partial\mathbf{h}}{\partial\mathbf{x}} {\frac{\partial f}{\partial\mathbf{x}}}^T = \frac{\partial\mathbf{h}}{\partial\mathbf{x}} {\frac{\partial\mathbf{h}}{\partial\mathbf{x}}}^T \boldsymbol{\lambda}.$$

and therefore by multiplying with the inverse of the non-singular matrix on the right side, we obtain:

$$- \left(\frac{\partial\mathbf{h}}{\partial\mathbf{x}} {\frac{\partial\mathbf{h}}{\partial\mathbf{x}}}^T\right)^{-1} \frac{\partial\mathbf{h}}{\partial\mathbf{x}} {\frac{\partial f}{\partial\mathbf{x}}}^T = \boldsymbol{\lambda}.$$

Their argument is that if this equation holds for an extremum at $$(\overline{\mathbf{x}}, \overline{\boldsymbol{\lambda}})$$, it approximately holds at $$\hat{\mathbf{x}} \approx \overline{\mathbf{x}}$$ and $$\hat{\boldsymbol{\lambda}} \approx \overline{\boldsymbol{\lambda}}$$, so we can insert $$\hat{\mathbf{x}}$$ and obtain $$\hat{\boldsymbol{\lambda}}$$.