Algorithms

Obtain bath fitting from pole fitting

In bath fitting, given $\mathrm{\Delta }(\mathrm{i}{\nu }_k)$ evaluated on $\{\mathrm{i}{\nu }_k{\}}_{k=1}^{N_w}$, we wish to find $V_j,E_j$ such that

$$\mathrm{\Delta }(\mathrm{i}{\nu }_k)=\sum _{j=1}^{N_b}\frac{V_j{V}_j^{†}}{\mathrm{i}{\nu }_k-E_j}.$$

This is achieved by the following strategy:

• Find pole fitting with semidefinite constraints:

  $$\begin{array}{}\text{(1)}& \mathrm{\Delta }(\mathrm{i}{\nu }_k)=\sum _{p=1}^{N_p}\frac{M_p}{\mathrm{i}{\nu }_k-{\lambda }_p},M_p\ge 0,\end{array}$$

  Here $M_p$ are $N_{\text{orb}}\times N_{\text{orb}}$ positive semidefinite matrices.

• Compute eigenvalue decomposition of $M_p$:

  $$\begin{array}{}\text{(2)}& M_p=\sum _{j=1}^{N_{\text{orb}}}{V}_j^{(p)}(V_j^{(p)}{)}^{†}.\end{array}$$

• Combining $\text{(1)}$ and $\text{(2)}$, we obtain the desired bath fitting:

$$\mathrm{\Delta }(\mathrm{i}{\nu }_k)=\sum _{pj}\frac{V_j^{(p)}(V_j^{(p)}{)}^{†}}{\mathrm{i}{\nu }_k-{\lambda }_p}.$$

Rational approximation via (modified) AAA algorithm

To find the poles ${\lambda }_p$ in $\text{(1)}$, we use the AAA algorithm, which is a rational approximation algorithm based on the Barycentric interpolant:

$$\begin{array}{}\text{(3)}& f(z)=\frac{p(z)}{q(z)}=\frac{\sum _{j=1}^m\frac{c_j{f}_j}{z-z_j}}{\sum _{j=1}^m\frac{c_j}{z-z_j}.}.\end{array}$$

The AAA algorithm is an iterative procedure. It selects the next support point in a greedy fashion. Suppose we have obtained an approximant $\tilde{f}$ from the $(k-1)$-th iteration, using support point $z_1,\cdots z_{k-1}$. At the $k$-th iteration, we do the following:

1. Select the next support point $z_k$ at which the previous approximant $\tilde{f}$ has the largest error.

2. Find $c_k$ in $\text{(3)}$ by solving the following linear square problem:

   $$\underset{\{c_k\}}{min}\sum _{z\ne z_1,\cdots z_k}{|f(z)q(z)-p(z)|}^2.\text{s.t.}\Vert c{\Vert }_2=1.$$

   This is a linear problem and amounts to solving a SVD problem. (See details in paper).

3. If the new approximant has reached desired accuracy, stop the iteration. Otherwise, repeat the above steps.

The poles of $f(z)$ are the zeros of $q(z)$, which can be found by solving the following generalized eigenvalue problem:

$$\begin{array}{r}\left(\begin{array}{ccccc}0& c_1& c_2& \cdots & c_m\\ 1& z_1& & & \\ 1& & z_2& & \\ ⋮& & & \ddots & \\ 1& & & & z_m\end{array}\right)=\lambda \left(\begin{array}{lllll}0& & & & \\ & 1& & & \\ & & 1& & \\ & & & \ddots & \\ & & & & 1\end{array}\right)\end{array}$$

For our application, we modify the AAA algorithm to deal with matrix-valued functions by replacing $f_j$ with matrices $F_j$.

Semidefinite programming

After obtaining ${\lambda }_p$, we need to find the weight matrices $M_p$ in $\text{(1)}$. We are solving the following problem:

$$\underset{\{M_p\}}{min}\sum _{k=1}^{N_w}{\Vert \mathrm{\Delta }(\mathrm{i}{\nu }_k)-\sum _{p=1}^{N_p}\frac{M_p}{\mathrm{i}{\nu }_k-{\lambda }_p}\Vert }_F^2,\text{s.t. }{M}_p\ge 0.$$

This is a linear problem with respect to $M_p$, and has semidefinite constraints, therefore could be solved efficiently via standard semidefinite programming (SDP) solvers.

Bi-level optimization

With ${\lambda }_p$ and $M_p$ obtained, we can further refine the poles and weights by solving the following bi-level optimization. Let us define the error function as

$$\text{Err}({\lambda }_p,M_p)=\sum _{k=1}^{N_w}{\Vert \mathrm{\Delta }(\mathrm{i}{\nu }_k)-\sum _{p=1}^{N_p}\frac{M_p}{\mathrm{i}{\nu }_k-{\lambda }_p}\Vert }_F^2.$$

Note that $\text{Err}$ is linear in $M_p$ but nonlinear in ${\lambda }_p$. As we have mentioned, optimization in $M_p$ is a SDP problem and therefore is robust, while optimization in ${\lambda }_p$ is a nonlinear problem and could be very challenging. This motivates us to define $\text{Err}({\lambda }_1,\cdots ,{\lambda }_{N_p})$ as a function of $\{{\lambda }_p\}$ only:

$$\text{Err}({\lambda }_1,\cdots ,{\lambda }_{N_p})=\underset{\{M_p\}}{min}\text{Err}({\lambda }_p,M_p)$$

The value of $\text{Err}({\lambda }_1,\cdots ,{\lambda }_{N_p})$ is obtained by solving a SDP problem. The gradient of $\text{Err}({\lambda }_1,\cdots ,{\lambda }_{N_p})$ could also be obtained analytically. (For details, see eq. 28 here.) And thus we could use a gradient-based optimization algorithm (L-BFGS) to minimize $\text{Err}({\lambda }_1,\cdots ,{\lambda }_{N_p})$ with respect to $\{{\lambda }_p\}$.

For performances, robustness and other details of this bi-level optimization framework, see again our original paper.