Algorithms

Obtain bath fitting from pole fitting

In bath fitting, given \(\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

\[\begin{equation} \Delta(\mathrm i\nu_k) = \sum_{j=1}^{N_b} \frac{V_jV_j^{\dagger}}{\mathrm i\nu_k - E_j}. \end{equation}\]

This is achieved by the following strategy:

  • Find pole fitting with semidefinite constraints:

    \[\begin{equation} \Delta(\mathrm i\nu_k) = \sum_{p=1}^{N_p} \frac{M_p}{\mathrm i\nu_k - \lambda_p}, M_p\geq 0, \tag{1} \label{polefit} \end{equation}\]

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

  • Compute eigenvalue decomposition of \(M_p\):

    \[M_p = \sum_{j=1}^{N_{\text{orb}}} V_{j}^{(p)} (V_{j}^{(p)})^{\dagger}. \tag{2} \label{eigdecomp}\]

  • Combining \(\eqref{polefit}\) and \(\eqref{eigdecomp}\), we obtain the desired bath fitting:

\[\Delta(\mathrm i\nu_k) = \sum_{pj} \frac{V_{j}^{(p)}(V_{j}^{(p)})^{\dagger}}{\mathrm i\nu_k - \lambda_p}.\]

Rational approximation via (modified) AAA algorithm

To find the poles \(\lambda_p\) in \(\eqref{polefit}\), we use the AAA algorithm, which is a rational approximation algorithm based on the Barycentric interpolant:

\[\begin{equation} f(z) = \frac{p(z)}{q(z)} = \frac{\sum_{j=1}^{m} \frac{c_jf_j}{z - z_j}}{\sum_{j=1}^{m} \frac{c_j}{z - z_j}.}. \tag{3} \label{bary} \end{equation}\]

The AAA algorithm is an iterative procedure. It selects the next support point in a greedy fashion. Suppose we have obtained an approximant \(\widetilde 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 \(\widetilde f\) has the largest error.

  2. Find \(c_k\) in \(\eqref{bary}\) by solving the following linear square problem:

    \[\begin{equation} \min_{\{c_k\}} \sum_{z\neq z_1,\cdots z_k} \left| f(z) q(z) - p(z) \right|^2. \quad \text{s.t.} \|c\|_2= 1. \end{equation}\]

    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{split}\begin{equation} \left(\begin{array}{ccccc} 0 & c_1 & c_2 & \cdots & c_m \\ 1 & z_1 & & & \\ 1 & & z_2 & & \\ \vdots & & & \ddots & \\ 1 & & & & z_m \end{array}\right)=\lambda\left(\begin{array}{lllll} 0 & & & & \\ & 1 & & & \\ & & 1 & & \\ & & & \ddots & \\ & & & & 1 \end{array}\right) \end{equation}\end{split}\]

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 \(\eqref{polefit}\). We are solving the following problem:

\[\begin{equation} \min_{\{M_p\}} \sum_{k=1}^{N_w} \left\| \Delta(\mathrm i\nu_k) - \sum_{p=1}^{N_p} \frac{M_p}{\mathrm i\nu_k - \lambda_p} \right\|_F^2, \quad \text{s.t. } M_p\geq 0. \end{equation}\]

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

\[\begin{equation} \text{Err}(\lambda_p, M_p) = \sum_{k=1}^{N_w} \left\| \Delta(\mathrm i\nu_k) - \sum_{p=1}^{N_p} \frac{M_p}{\mathrm i\nu_k - \lambda_p} \right\|_F^2. \end{equation}\]

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:

\[\begin{equation} \text{Err}(\lambda_1,\cdots, \lambda_{N_p}) = \min_{\{M_p\}}\text{Err}(\lambda_p, M_p) \end{equation}\]

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.