TITLE: When does a Hermitian operator have real matrix elements?
QUESTION [3 upvotes]: I will use braket-notation, but my question is not specific to quantum mechanics. Instead, I would be interested in a general answer for operators in some Hilbert space. Let $H$ be a Hermitian operator with eigenstates $|i\rangle$, so that $H |i\rangle = E_i |i\rangle$, where some eigenvalues may possibly be degenerate. Now consider another Hermitian operator $A$. This operator can be represented as a matrix in the basis $\{|i\rangle\}$ of the eigenvectors of $H$, with elements
$$
A_{ij} = \langle i|A|j \rangle
$$
Hermiticity of $A$ then requires $A_{ji} = A_{ij}^*$. In general, however, these matrix elements may be complex. My question is the following: Is it possible to formulate a condition on $A$, probably in relation to $H$, such that the matrix elements of $A$ in the basis formed by the eigenvectors of $H$ are real, $A_{ji} = A_{ij}$, if this condition is satisfied?
I think in some situations a simple multiplication by a phase factor may be sufficient. Assuming the $A_{ij}$ are complex, one may write
$$
A_{ij} = |A_{ij}| \, e^{i \phi_{ij}}
$$
Now consider transforming to new basis vectors given by $|i'\rangle = e^{i \nu_i} |i\rangle$. These are still eigenstates of $H$ and the matrix elements of $A$ in this new basis are given by
$$
A_{ij}' = |A_{ij}| \, e^{i (\phi_{ij} + \nu_j - \nu_i)}
$$
So if there is a solution for the $\nu_i$ to the set of $n^2$ equations (where $n$ is number of eigenstates of $H$, so dimension of Hilbert space) given by
$$
\phi_{ij} + \nu_j - \nu_i = 0,
$$
then the operator $A$ can be represented by a real matrix in that basis. I believe that such a solution exists in the case when the phases satisfy the relation $\phi_{ij} + \phi_{jk} = \phi_{ik}$. However, I don't think the phases of a general operator $A$ necessarily satisfy this condition. If they do not, there may not be a solution to the system of equations since there are $n^2$ constraints, but only $n$ variables $\nu_i$ to solve for.
Is there a general relation between $A$ and $H$ that leads to a representation of $A$ in terms of a matrix with real elements?

REPLY [0 votes]: I don't think you will get $n^2$ unique equations, since only off-diagonal elements can be complex. Further you need to only consider the upper off-diagonal ones since the bottom off-diagonals are just their complex conjugates. That leaves you with a total set of $\frac{n(n-1)}{2}$ phases $\phi_{ij}$. Rescaling the $n$ basis vectors by an arbitrary phase also gives you $\frac{n(n-1)}{2}$ possibilities for the differences $\nu_i - \nu_j$.
For example, take the case $n=3$. We have the phases $\phi_{12}, \phi_{13}, \phi_{23}$ for the upper off-diagonal elements and minus those for the lower off-diagonal elements. Rescaling the basis elements, we have
$$
\phi_{12} = \nu_1 - \nu_2,\\
\phi_{13} = \nu_1 - \nu_3, \\
\phi_{23} = \nu_2 - \nu_3.
$$
EDIT:
As user @Dvij D.C. pointed out, this case works since $\frac{n(n-1)}{2}=n$ has a solution when $n=3$, i.e. the number of independent variables is equal to the number of equations for the differences. In higher dimensions, the number of equations will be more than the number of independent variables and so in general this will not be possible unless the phases are chosen in a specific way.


TITLE: Does a free electron, one that's not either in an atom or a wire, have an associated wave-function?
QUESTION [2 upvotes]: Would a free electron, one that's not either in an atom or moving through a wire, but moving through empty space on its own, have an associated wave-function? Or, is an electron described as a wave-function only when inside an atom?
And if the type of free electron I'm describing does have a wave-function, how can it when the electron wouldn't be in an orbital within an atom?

REPLY [1 votes]: There is no need of a potential for the Schrödinger equation to have a solution, namely
$$
i\hbar \frac{\partial}{\partial t} |\psi\rangle = H |\psi\rangle
$$
does possess solutions even when the Hamiltonian contains no potential. Questions of normalisability may arise, but that is another point (and can anyway be solved by expanding in Fourier terms).

And if the type of free electron I'm describing does have a wave-function,   how can it when the electron wouldn't be in an orbital within an atom?

You seem to have a little odd idea of what a wave function is, probably due to some odd-written books in chemistry which associate the wave function to some sort of orbital cloud. That is by no means what the wave function is. Rather it is the representation on the position basis  $\langle x|\psi\rangle$ of the solution of the Schrödinger equation above.


TITLE: Linear Programming: "at most k out of n variables nonzero" constraint
QUESTION [4 upvotes]: I have a Mixed-Integer Program that contains (among other things) $n$ variables $v_1, \dots v_n$ (continuous or integer doesn't matter, in $[0, M)$ for some $M$).
I want to formulate the constraint that at most $k < n$ of these variables may become nonzero. The obvious solution would be to introduce a set of binary indicator variables $w_1, \dots w_n$, a big-M constraint per pair and a sum constraint on the indicators:
$$ w_i \cdot M \geq v_i$$
$$ \sum w_i \leq k$$
However, that's a lot of big-M's and probably wreaks havoc on the relaxed solution. I'm wondering if there are standard techniques for coping with a setting like this.
What I've thought of: Introduce a SOS-Type 1 - Constraint per $v_i$ / $w_i$ pair. This would get rid of the big-M and force $w_i$ to become zero if $v_i$ is nonzero. (Obviously, the second constraint above must be changed to $\sum (1 - w_i) \leq k$). Is that the way to go? It feels like I'm using SOS-constraints for something that they weren't meant to (it's not an ordered set…)
Any help is appreciated! Thanks, Lukas

REPLY [3 votes]: Internally, a solver is probably going to treat each of your SOS1 constraints a lot like a binary variable -- basically, branch on $w_i$ and, in the $w_i = 1$ branch, force $v_i = 0$. The good news is that you avoid $M$; the bad news is that I don't think the solver has much information with which to assess which SOS1 constraint to branch on at a given node.
If the variables have natural (and not horribly large) upper bounds, the big $M$ approach may not to be too bad. If you don't know those bounds, you might relax integrality conditions and maximize each variable in turn to get bounds.
Another possibility is combinatorial Benders decomposition (Codato and Fischetti, Operations Research 54 (2006), 756-766). Rather than use big $M$ constraints, you have a master problem containing the $w_i$ and one continuous variable (a surrogate for the objective contribution of the $v_i$), and a subproblem that optimizes the original objective function over the $v_i$ subject to the original constraints. For a given solution to the master problem, you simply omit $v_i$ from the subproblem if $w_i = 0$. If the subproblem is infeasible, you add a "feasibility cut" to the master. If the subproblem has an optimal solution but its objective value is worse than the value of the surrogate variable in the master solution, you add a constraint ("optimality cut") to the master problem tying the value of the surrogate to the values of the $w_i$.
I would only use combinatorial Benders if the $M$ formulation was either seriously slow or had numerical stability problems (meaning, in either case, that the "big $M$" coefficients really were big). Also, I'd be hesitant to do it if some or all of the $v_i$ are integer valued, since (a) you'd be repeatedly solving a MIP subproblem and (b) the optimality and feasibility cuts you got might be weaker than usual.


TITLE: Binomial expansion with 2n power
QUESTION [2 upvotes]: I'm trying to expand $$(z+z^{-1})^{2n}$$.
I've written the first 3 terms out but I can't figure out a way to work out the last 3 terms.
I'm trying to produce an expression for $$\cos^{2n}x$$ using the complex numbers.

REPLY [1 votes]: Just apply the binomial formula
$$
 (z+z^{-1})^{2n}=\sum_{k=0}^{2n}\binom{2n}kz^k(z^{-1})^{2n-k}
 =\sum_{k=0}^{2n}\binom{2n}kz^{2(k-n)}
 =\sum_{k=-n}^n\binom{2n}{n+k}z^{2k}.
$$
Now putting $\def\i{\mathbf i}z=\exp(\i x)$ and dividing by $2^{2n}$ you get
$$
  \cos(x)^{2n}= 2^{-2n}\sum_{k=-n}^n\binom{2n}{n+k}\exp(2\i kx)
 = 2^{-2n}\sum_{k=-n}^n\binom{2n}{n+k}\cos(2kx).
$$


TITLE: Expected value of minimum distance
QUESTION [3 upvotes]: I'm trying to estimate average complexity of some algorithm and the problem I'm currently stuck with can be stated very simply.
Let $a_i$ be a list of $m$ uniformly distributed integers not greater than some $n$: $a_i \in \mathbb{N} \land a_i \leq n$ for $i={1..m}$.
Then I pick some $j$ randomly and calculate minimum distance between $a_j$ and other $a_i$'s: $x = \min_{j \neq i} |a_j - a_i|$.
The question is how to calculate expected value of $x$?

REPLY [1 votes]: Answer: $\left(\frac{1}{2}+o(1)\right) \frac{n}{m} + \delta$, where $|\delta| < 1$ and $o(1)$ tends to $0$ as $m\to \infty$.
First of all, note that since all $a_i$ are independently and identically distributed,
$$\mathbb{E}[x| j =1] = \mathbb{E}[x| j = 2] = \dots = \mathbb{E}[x| j =m],$$
so it doesn't matter whether we choose $j$ randomly or just let $j=1$. Let's assume that $j=1$.
It's easier to work with continuous random variables. So let's assume that we choose independent random numbers $y_1, \dots, y_m \in_U (0,1)$ and then let $x_i = \lceil ny_i\rceil$. Note that
$$n|y_i - y_j| - 1 < |x_i - x_j| < n|y_i - y_j| + 1.$$
Therefore, if we are only interested in an approximate answer, it suffices to find the expectation $\mathbb{E}[|y_i - y_j|]$. We reduced our problem to the following problem.

We are given $m$ independent random variables $y_1, \dots, y_m$ uniformly distributed on $(0,1)$. Compute the expectation of $\xi = \min_{i>1} |y_1-y_i|$.

Here is an informal solution to this problem.
\begin{align*}\mathbb{E}[\xi] &= \int_0^1 \Pr[\xi > t] dt =  \int_0^1 \Pr[|y_1 - y_i| > t \text{ for every } i > 2]\, dt \\ &= \int_0^1 \mathbb{E}[\Pr[|y_1 - y_2| > t]^{m-1}|y_1] \,dt.
\end{align*}
Since there are $m$ numbers $\{y_i\}$ on $(0,1)$, we expect that the distance between consecutive numbers will be about $1/m$; it's very unlikely that the distance is more than $\tau = \log m /m$. Thus $\xi < \tau$ with very high probability, also $y_1 \in(\tau, 1 - \tau)$ w.h.p. Let us condition on $y_1 = y^*$ where $y^* \in (\tau, 1 - \tau)$. Then we have
\begin{align*}
\int_0^1 \Pr[|y_1 - y_2| > t | y_1 = y^*]^{m-1} dt &\approx \int_0^{\tau} \Pr[|y_2 - y^*| > t]^{m-1} dt \\
&= \int_0^{\tau} \Pr[y_2\notin [y^* -t, y^*+t]]^{m-1} dt\\
&=\int_0^{\tau} (1 - 2t)^{m-1} dt.
\end{align*}
Since $y_1 \in(\tau, 1 - \tau)$ w.h.p, we have
$$\mathbb{E}[\xi] \approx \int_0^{\tau} (1 - 2t)^{m-1} dt \approx \frac{1}{2m}.
$$
Finally, $\mathbb{E}[\xi] \in (n\mathbb{E}[\xi] - 1, n\mathbb{E}[\xi] +1)$. Plugging in the expression for $\mathbb{E}[\xi]$, we get
$$\mathbb{E}[\xi] \in \left(\left(\frac{1}{2}+o(1)\right) \frac{n}{m} -1 , \left(\frac{1}{2}+o(1)\right) \frac{n}{m}  +1\right).$$
Note that this bound is not particularly useful when $m\approx n$.


TITLE: Question regarding congruence equation of the form $ax^2 + bx + c \equiv 0 \pmod{m}$
QUESTION [3 upvotes]: Problem
  Prove that the congruence equation $ax^2 + bx + c \equiv 0 \pmod{m}$, with $(2a, m) = 1$ is equivalent to a congruence of the form $x^2 \equiv r \pmod{m}$.

I really have no idea where to start. Is there a theorem which relates to this problem? What do we have to show? A hint would be sufficient.
Thanks,

REPLY [2 votes]: Multiply it by the unit $\rm\:4a\:$ yielding the equivalent equation $\rm\ (2\:a\:x+b)^2\: =\ b^2 - 4\:a\:c\ \ (mod\ m)$

REPLY [1 votes]: Remember the quadratic formula to solve $ax^2 + bx + c = 0$: $x = \tfrac{1}{2a}(-b \pm \sqrt{b^2 - 4ac})$.
This works for the congruence too if we just interpret it right.
Suppose $r^2 \equiv b^2 - 4ac$ then it's easy to check that $x \equiv (2a)^{-1} (-b \pm r)$ solves the quadratic. The condition 2a coprime is m is exactly what's needed to invert 2a.
if x solves the quadratic congruence then multiply by 2a and add b to get the other direction.


text:
\section{Anytime pMHE Algorithm}
\label{sec: iteration scheme}
In this section, we propose a novel pMHE iteration scheme in which, at each time instant $k$, problem (\ref{condensed}) is approximately solved by executing  a fixed number $\mathrm{it}(k) \in \mathbb{N}_+$ of optimization algorithm iterations. In more details, at each time $k$, a suitable warm start $\mathbf{\hat{z}}^{0}_{k}$ is generated from a stabilizing a priori estimate $\mathbf{\bar{z}}_k$ and an iterative optimization update is  carried out, from which the sequence   $\big\lbrace  \mathbf{\hat{z}}^{i}_{k}   \big\rbrace $ with $i=1,\cdots,\mathrm{it}(k)$ is obtained. The steps of the scheme are given in Algorithm~\ref{alg: eta changes} and illustrated in Figure \ref{Fig: iteration scheme}.
\begin{figure}[ht]
  \centering
\input{figures/iteration_scheme_sketch}
\caption{Illustration of the steps of the pMHE iteration scheme at the time instants $k-1$ and $k$, with the corresponding $\mathrm{it}(k-~1)=4$  and $\mathrm{it}(k)=3$ optimization iterations. }
\label{Fig: iteration scheme}
\end{figure}
~\\Before we explain the proposed  algorithm in more detail, and for the sake of clarity, let us first introduce  some notations. The index $k$ denotes the time instant in which we receive a new measurement and $\mathrm{it}(k)$ is the number of iterations of the optimization algorithm between time instants $k$ and $k+1$. Moreover, we introduce
\begin{align}
\label{notation}
\hspace{-0.4cm}\mathbf{\hat{z}}^{ i}_{k} \coloneqq \begin{bmatrix}
				 \hat{x}^{ i}_{k-N} \\[0.2em]  \hat{w}_{k-N}^{ i} \\\vdots  \\[0.2em]  \hat{w}_{k-1}^{ i}    \end{bmatrix}, \enskip
 \mathbf{\bar{z}}_k \coloneqq \begin{bmatrix}
 \bar{x}_{k-N}\\[0.2em]  \bar{w}_{k-N}\\\vdots  \\[0.2em]  \bar{w}_{k-1}
 \end{bmatrix}, \enskip
     \mathbf{z}_k \coloneqq \begin{bmatrix}
				 x_{k-N}  \\[0.21em] 0 \\\vdots  \\[0.21em]  0  \end{bmatrix}.
 \end{align}
  With $\mathbf{\hat{z}}^{ i}_{k} $, we denote the $i$-th iterate of the optimization algorithm at time $k$.
With $ \mathbf{\bar{z}}_k $, we refer to the a priori estimate at time $k$ and with   $ \mathbf{z}_k$ 	 to the true state $x_{k-N}$ with true model residual sequence    $  \lbrace  0  , \cdots , 0  \rbrace $.

Upon arrival of a new measurement at time $k$, the optimization algorithm is initialized based on the a priori estimate~$\mathbf{\bar{z}}_k$.  In particular, we compute the warm start $\mathbf{\hat{z}}^{0}_{k}$ as the Bregman projection  of $\mathbf{\bar{z}}_k$ onto the constraint set $ \mathcal{S}_k$ as formulated in line~3 of Algorithm \ref{alg: eta changes}.
Then, a fixed number $\mathrm{it}(k)$ of optimization iterations is performed via (\ref{optAdap}), generating $\big\lbrace  \mathbf{\hat{z}}^{1}_{k}, \cdots,\mathbf{\hat{z}}^{\mathrm{it}(k)}_{k} \big\rbrace $. Here, $\eta_k^i >0$ denotes the step size employed at the $i$-th iteration at time $k$. From this sequence of iterates, an arbitrary iterate  $\mathbf{\hat{z}}^{\mathrm{j}(k)}_{k} $, $\mathrm{j}(k) \in \lbrace 0,\cdots,\mathrm{it}(k)\rbrace$ can be chosen, based on which the state estimate $\hat{x}_k$ is obtained using
\begin{align}
		\label{StateEst}
\hat{x}_k =A^N\, \hat{x}^{\mathrm{j}(k)}_{k-N}+ \!\!\sum_{j=k-N} ^ {k-1}\! \! A^{k-1-j}\, \left( B\,  u_j+ \hat{w}^{\mathrm{j}(k)}_{j} \right)
\end{align}
for $k>N$ (see Remark \ref{rem: indices} for the case where $0< k \leq N$). Moreover, the a priori estimate $ \mathbf{\bar{z}}_{k+1} $ for the next time instant is computed through the operator $\Phi_k : \mathbb{R}^{(N+1)n} \rightarrow \mathbb{R}^{(N+1)n}$ which will be defined in~\eqref{Adap: warm}.  As mentioned above, the basic idea of the  pMHE framework is to use the Bregman distance~$D_\psi$ as a proximity measure to a  stabilizing a priori estimate in order to inherit its  stability properties. Since the Luenberger observer  appears as a simple candidate for constructing the a priori estimates,  we require that the operator $\Phi_k$ incorporates its dynamics as follows:
			    \begin{align}
			     \label{Adap: warm}
\Phi_k \left( \mathbf{\hat{z}}^{ \mathrm{j}(k)}_{k} \right) \! \coloneqq  \!  \begin{bmatrix}
   \! A\,  \hat{x}_{k-N}^{ \mathrm{j}(k)} +\!B\, u_{k-N}+\!L \left(y_{k-N} - C  \hat{x}_{k-N}^{ \mathrm{j}(k)} \right)     \\\mathbf{0}
     \end{bmatrix}\!,
       \end{align}
where $\mathbf{0}\in \mathbb{R}^{Nn}$. Here, the observer gain $L$ is chosen such that all the eigenvalues of $A-LC$ are strictly within the unit circle.
 \begin{algorithm} [t]
\caption{Anytime pMHE}
\label{alg: eta changes}
\begin{algorithmic}[1]
\State  \textbf{Initialize:}   Choose  $ \hat{x}_0  $ and set $ \bar{\mathbf{z}}_0  =\hat{x}_0$
  \For {$k=1,2, \cdots$  }
    \State  $ \mathbf{\hat{z}}_k^0= \argmin\limits_{\quad\enskip\mathbf{z}\in  \mathcal{S}_k} \quad  D_\psi (\mathbf{z},\mathbf{\bar{z}}_k)$  warm start
      \For {$i= 0,\dots, \mathrm{it}(k)-1$ } optimizer update
      \begin{align}
  \label{optAdap}
\hspace*{-0.2cm}  \mathbf{\hat{z}}^{i+1}_{k}=\argmin_{\quad\enskip\mathbf{z}\in  \mathcal{S}_k}\left\lbrace  \eta_k^i\,    \nabla f_k\left(  \mathbf{\hat{z}}^{i}_{k} \right)^\top \mathbf{z}+ D_\psi (\mathbf{z},  \mathbf{\hat{z}}^{i}_{k}) \right\rbrace
   \end{align}
		\EndFor
		\State for some $\mathrm{j}(k) \in \lbrace 0,\cdots,\mathrm{it}(k)\rbrace$ obtain $\hat{x}_k$ from (\ref{StateEst})
		\State  $\mathbf{\bar{z}}_{k+1} =\Phi_k \left( \mathbf{\hat{z}}^{\mathrm{j}(k)}_{k} \right)$
		\EndFor
\end{algorithmic}
\end{algorithm}
In the following, we compare Algorithm~\ref{alg: eta changes} with our earlier formulation of pMHE, given in  Algorithm~\ref{alg: optimal pMHE}.
 Observe that, while a solution of the optimization problem  (\ref{condensedpMHE}) is computed,  step~4 in Algorithm \ref{alg: eta changes} employs the so-called  mirror descent algorithm~\cite{beck2003mirror} that iterates~\eqref{optAdap}  until  a given number of iterations $\mathrm{it}(k)$ is achieved. For $D_\psi \left( \mathbf{z}_1, \mathbf{z}_2 \right) = \frac{1}{2} \left\Vert \mathbf{z}_1 -  \mathbf{z}_2 \right\Vert^2$ and $ \mathcal{S}_k=\mathbb{R}^{(N+1)n}$, the optimizer update step (\ref{optAdap}) corresponds to an iteration step of the classical gradient descent algorithm and hence step 4
 can be executed very quickly.
In the constrained case, (\ref{optAdap}) can be regarded as a generalization of the projected gradient algorithm \cite{beck2003mirror}.  For this reason, we can view Algorithm  \ref{alg: eta changes} as a real-time version of the pMHE scheme given  in  Algorithm~\ref{alg: optimal pMHE}. \rev{Note also that choosing the so-called Kullback-Leiber divergence  as Bregman distance yields the efficient entropic descent algorithm if the constraint set is given by the unit simplex  \cite{beck2003mirror}.} An appealing feature of  Algorithm \ref{alg: eta changes} is that, depending on the available computation time between two subsequent time instants $k$ and $k+1$, the user can specify a maximum number of iterations $\mathrm{it}(k)$ after which   the optimization algorithm at time $k$ has to return a solution. \\
Another key difference between the two algorithms is that Algorithm \ref{alg: optimal pMHE}  is biased by the stabilizing a priori estimate $\bar{\mathbf{z}}_k$, while this bias  is fading away  in  Algorithm \ref{alg: eta changes}. In other words, the a priori estimate constructed based on the Luenberger observer (\ref{Adap: warm}) has less impact at each optimization iteration, which improves the performance of the pMHE iteration scheme.
This is due to the fact that $\bar{\mathbf{z}}_k$
might  degenerate  performance in Algorithm \ref{alg: optimal pMHE}, since
the  solution lies in proximity to the a priori estimate.
In Algorithm \ref{alg: eta changes}, however, the Luenberger observer enters
only in the warm start. From this point of view,
it is quite surprising that, even though the effect of this stabilizing ingredient is fading away, stability is provably preserved, as we will show in the subsequent section. Thus, this "implicit stabilizing regularization"
approach of the a priori estimate
is in contrast to the explicit stabilizing regularization proposed in \cite{sui2014linear,gharbi2018proximity} (see also \cite{belabbas2020implicit}).
Moreover, the proposed MHE algorithm possesses the anytime property.
  The \emph{anytime} property refers to the fact that the algorithm will yield stable estimation errors after any number of  optimization algorithm iterations. This is similar in spirit to anytime model predictive control (MPC) algorithms, which compute stabilizing  control inputs after any  optimization iteration \cite{feller2017stabilizing,feller2018sparsity}.
\begin{rem} \label{rem: indices} For $k \leq N$, we can employ the steps of Algorithm \ref{alg: eta changes} by setting all the negative indices to zero. More specifically,
$  \mathbf{\hat{z}}^{ i}_{k} =\begin{bmatrix} (\hat{x}^{ i}_{0})^\top &  (\hat{w}_{0}^{ i})^\top & \hdots  & (\hat{w}_{k-1}^{ i})^\top    \end{bmatrix}^\top$ and  $\mathbf{\bar{z}}_{k}  =\begin{bmatrix} \hat{x}_{0}^\top &  0 & \hdots  & 0 \end{bmatrix}^\top$ where $\mathbf{\hat{z}}^{ i}_{k},\mathbf{\bar{z}}_{k}   \in \mathbf{R}^{(k+1)n} $. In order to compute the state estimate, \eqref{StateEst} has to be explicitly modified to
 \begin{align}
\hat{x}_k =A^k\, \hat{x}^{ \mathrm{j}(k)}_{0}+  \sum_{j=0} ^ {k-1}  A^{k-1-j}\,\left( B\,   u_j + \hat{w}^{ \mathrm{j}(k)}_{j} \right).
\end{align}
 \end{rem}
We impose the following assumptions.
\begin{ass}[Properties of $ \mathcal{S}_k$] \label{ass: S closed}
 The  set $ \mathcal{S}_k$ of constraints is closed and convex with nonempty interior.
\end{ass}
  \begin{ass}[Convexity of $f_k$] \label{ass: f convex}
The sum of stage costs  $f_k$ is continuously differentiable,  convex for all $k> 0$, and achieves its minimum at $\mathbf{z}_k$.
   \end{ass}
     \begin{ass}[Strong smoothness of $f_k$] \label{ass: f strngSmooth}
The sum of stage costs  $f_k$ is strongly smooth with constant $L_f>0 $:
\begin{align}
 f_k(\mathbf{z}_2) \leq  f_k(\mathbf{z}_1)  +\nabla f_k(\mathbf{z}_1)^\top   ( \mathbf{z}_2 -\mathbf{z}_1)  + \frac{L_f}{2} \Vert  \mathbf{z}_1 - \mathbf{z}_2  \Vert^2
\end{align}
for all $\mathbf{z}_1, \mathbf{z}_2 \in \mathbb{R}^{(N+1)n}  $ and $k>0$.
   \end{ass}
 \begin{ass} [Strong convexity and smoothness of $D_\psi$] \label{ass: Bregman}  The function $\psi$ is continuously differentiable, strongly convex  with constant $\sigma>0$ and strongly smooth with constant  $\gamma>0$, which implies the following for the Bregman distance
 \begin{align}
 \frac{\sigma}{2 } \Vert  \mathbf{z}_1 - \mathbf{z}_2  \Vert^2  \leq   D_\psi (\mathbf{z}_1, \mathbf{z}_2 ) \leq  \frac{\gamma}{2} \Vert  \mathbf{z}_1 - \mathbf{z}_2  \Vert^2
 \end{align}
 for all $\mathbf{z}_1, \mathbf{z}_2 \in \mathbb{R}^{(N+1)n}  $.
   \end{ass}
In Assumption \ref{ass: f convex}, we can ensure that $f_k$ achieves its minimum at $\mathbf{z}_k$ by designing the stage costs $r(\cdot)$ and $q(\cdot)$ such that their corresponding minimum is achieved at zero. \rev{Requiring strong smoothness of $f_k$ in Assumption \ref{ass: f strngSmooth} is	rather customary in the analysis of first-order optimization algorithms \cite{mokhtari2016online} and will prove central in the subsequent theoretical studies of the pMHE iteration scheme. Notice that Assumption \ref{ass: Bregman} restricts the  class of employed  Bregman distances to functions which can be quadratically lower and upper bounded. Obviously, this includes the important special case of quadratic distances which are widely used as prior weighting in the MHE literature in order to ensure stability of the estimation error. In addition,we can use any Bregman distance $D_{\psi}$ constructed based on the function  $\psi(\mathbf{z})= \frac{1}{2}\Vert \mathbf{z} \Vert_{P}^2 + B(\mathbf{z})$, where $B:\mathbb{R}^{(N+1)n} \rightarrow \mathbb{R}$ is a convex and strongly smooth function, i.e., a convex function whose gradient is Lipschitz continuous. }

config: arxiv


TITLE: proof of l'Hôpital's rule that minimizes special-casing
QUESTION [3 upvotes]: A simple form of l'Hôpital's rule looks like this: If $u$ and $v$ are functions with $u(0)=0$ and $v(0)=0$, the derivatives $\dot{v}(0)$ and $\dot{v}(0)$ are defined,
and the derivative $\dot{v}(0)\ne 0$,
then
\begin{align*}
  \lim_{x\rightarrow 0} \frac{u}{v} &= \frac{\dot{u}(0)}{\dot{v}(0)} \qquad .
\end{align*}
To me, the clearest way to arrive at this result uses a little nonstandard analysis: Since $u(0)=0$, and the derivative $d u/d x$ is defined at $0$, $u(d x)=d u$ is infinitesimal, and likewise for $v$. By the definition
of the limit, the limit is the standard part of
\begin{equation*}
  \frac{u}{v} = \frac{d u}{d v} = \frac{d u/d x}{d v/d x} \qquad ,
\end{equation*}
where by assumption the numerator and denominator are both defined (and finite, because the derivative
is defined in terms of the standard part). The standard part of a quotient like $p/q$ equals the quotient of the
standard parts, provided that both $p$ and $q$ are finite (which we've established), and $q \ne 0$ (which is
true by assumption). But the standard part of $d u/d x$ is the definition of the derivative $\dot{u}$, and
likewise for  $d v/d x$, so this establishes the result.
The generalizations to $x\rightarrow a$, where $a\ne 0$, and $x\rightarrow \infty$ are pretty trivial with the changes of variable $x\rightarrow x-a$ and $x\rightarrow 1/x$.
But there are a bunch of other cases of l'Hôpital's that seem to me to involve toxic doses of case-splitting. There are cases where you have to differentiate more than once, and cases where the indeterminate form is $\infty/\infty$ rather than $0/0$.
Is it possible to treat all of this in a unified way, possibly using ideas from projective geometry or inversions with respect to a circle in the complex plane?

REPLY [3 votes]: Last week I covered L'Hôpital's Rule in my Spivak Calculus course.  I ended up preferring the proof given in Rudin's Principles to Spivak's treatment, because (i) he does compile several different cases in a rather efficient way and (ii) he proves a stronger result in the case where $\lim_{x \rightarrow a} g(x) = \infty$, namely that there is no hypothesis needed on the limiting behavior of the numerator $f$.
I wrote this up in $\S$ 1 of these notes.  However I must confess that I found the proof itself not very interesting and ended up not covering it in class.
I must also confess that I don't really understand the OP's sketch proof.  In particular I do not recognize the hypothesis $g'(a) \neq 0$ that the OP is assuming: who says that $f$ and $g$ are even defined at $a$?  I am also slightly skeptical that NSA really helps here to give a shorter proof, but I would be very interested to be proven wrong about this.
Added: If you are willing to assume that $f$ and $g$ are defined and differentiable at $a$, that $f(a) = g(a) = 0$ and that $g'(a) \neq 0$, the proof becomes almost trivial:
$\lim_{x \rightarrow a} \frac{f(x)}{g(x)} = \lim_{x \rightarrow a} \frac{f(x) - f(a)}{g(x)-g(a)} = \lim_{x \rightarrow a} \frac{ \frac{f(x)-f(a)}{x-a}}{\frac{g(x)-g(a)}{x-a}} = \frac{ \lim_{x \rightarrow a} \frac{f(x) - f(a)}{x-a}}{\lim_{x \rightarrow a} \frac{g(x)-g(a)}{x-a}} = \frac{f'(a)}{g'(a)}$.
(However this version is inadequate for many of the standard applications of freshman calculus.)  So I presume we're talking about a stronger version than this?