pymor.algorithms.loewner

Module Contents

pymor.algorithms.loewner.complete_conjugate_pairs(nodes, *data)

Complete nodes and associated data with complex conjugate pairs.

For each node whose complex conjugate is missing, append its conjugate and the elementwise conjugates of the corresponding entries in each array in data. Existing entries retain their order; new entries are appended in input-node order. For example, complete samples and weights together with:

nodes, samples, weights = complete_conjugate_pairs(nodes, samples, weights)

Note

Nodes are compared using exact equality. Already supplied conjugate data are not checked for consistency. Real nodes are not duplicated, and input arrays are not modified.

Parameters:

• nodes – Nonempty one-dimensional NumPy array of sampling nodes.

• data – Zero or more NumPy arrays passed as positional arguments, e.g. samples, derivatives, quadrature weights or tangential directions. Each array must have shape (len(nodes), ...); trailing dimensions are preserved.

Returns:

• completed_nodes – One-dimensional NumPy array containing the original and appended nodes.

• completed_data – Completed NumPy arrays in input order, returned as separate entries in (completed_nodes, *completed_data), not a nested tuple. With no data arguments, the return value is (completed_nodes,).

pymor.algorithms.loewner.loewner_matrices(left_nodes, right_nodes, left_terms, right_terms, derivative_terms=None)

Construct Loewner and shifted Loewner matrices from pairwise terms.

For distinct left and right nodes, the entries are

\[\mathbb{L}_{ij} = \frac{A_{ij} - B_{ij}}{\mu_i - \lambda_j}, \qquad (\mathbb{L}_s)_{ij} = \frac{\mu_i A_{ij} - \lambda_j B_{ij}}{\mu_i - \lambda_j}.\]

The shifted matrix is computed as \((\mathbb{L}_s)_{ij} = \mu_i\mathbb{L}_{ij} + B_{ij}\). At coincident nodes, this gives the Hermite value \(\mu_i D_{ij} + B_{ij}\), where \(D_{ij}\) is the supplied derivative term.

Parameters:

• left_nodes – Nonempty one-dimensional NumPy array of left interpolation nodes \(\mu_i\), of shape (n_left,).

• right_nodes – Nonempty one-dimensional NumPy array of right interpolation nodes \(\lambda_j\), of shape (n_right,).

• left_terms – NumPy array containing \(A_{ij}\). Together with right_terms, must broadcast to shape (n_left, n_right, ...).

• right_terms – NumPy array containing \(B_{ij}\), broadcastable with left_terms as above.

• derivative_terms – Optional NumPy array containing \(D_{ij}\), broadcastable to the pairwise term shape. Required at coincident nodes, where the left and right terms must agree. See loewner_matrix for the derivative and coincidence conventions.

Returns:

• L – Loewner matrix as a NumPy array of shape (n_left, n_right, ...).

• Ls – Shifted Loewner matrix as a NumPy array of the same shape as L.

pymor.algorithms.loewner.loewner_matrix(left_nodes, right_nodes, left_terms, right_terms, derivative_terms=None)

Construct a Loewner matrix from pairwise terms.

The returned array contains the divided differences

\[\mathbb{L}_{ij} = \frac{A_{ij} - B_{ij}}{\mu_i - \lambda_j}.\]

Scalar, matrix-valued and tangentially projected terms are supported through broadcasting. Use loewner_quadruple to assemble pairwise terms directly from transfer function data.

Note

Coincident nodes are detected using exact equality. Nearly coincident nodes are treated by divided differences and may suffer from cancellation.

Parameters:

• left_nodes – Nonempty one-dimensional NumPy array of left interpolation nodes \(\mu_i\), of shape (n_left,).

• right_nodes – Nonempty one-dimensional NumPy array of right interpolation nodes \(\lambda_j\), of shape (n_right,).

• left_terms – NumPy array containing \(A_{ij}\). Together with right_terms, must broadcast to shape (n_left, n_right, ...). For scalar function samples, a column of left values of shape (n_left, 1) can be supplied.

• right_terms – NumPy array containing \(B_{ij}\), broadcastable with left_terms as above. For scalar function samples, a row of right values of shape (1, n_right) can be supplied.

• derivative_terms – Optional NumPy array of derivative terms with respect to the complex argument, broadcastable to the pairwise term shape. Required when a left and right node coincide. At these entries, the left and right terms must agree and the supplied derivative replaces the divided difference. Other derivative entries are ignored.

Returns:

L – NumPy array of shape (n_left, n_right, ...) containing the Loewner matrix. Trailing dimensions of the broadcast terms are retained, not flattened into blocks.

pymor.algorithms.loewner.loewner_matrix_nd(sampling_values, samples, interpolation_indices)

Construct a higher-dimensional Loewner matrix for Cartesian scalar data.

Assemble the Loewner matrix used in the parametric AAA algorithm of [CRBG23]. The interpolation nodes form a Cartesian subgrid of the sampling grid. For one variable, this reduces to loewner_matrix with interpolation nodes on the right and all remaining nodes on the left.

Parameters:

• sampling_values – Nonempty list or tuple of nonempty one-dimensional NumPy arrays, one per variable. The first variable is typically the complex frequency; subsequent variables are parameters. Interpolation and non-interpolation node values must not coincide in any variable.

• samples – NumPy array of scalar samples with shape tuple(map(len, sampling_values)). samples[i, j, ...] corresponds to the nodes (sampling_values[0][i], sampling_values[1][j], ...). Matrix-valued samples must be projected to scalar data before calling this function.

• interpolation_indices – List or tuple of nonempty one-dimensional integer NumPy arrays, one per variable. Each array selects interpolation nodes from the corresponding sampling_values array and must contain distinct, valid indices. For one variable, at least one sampling node must remain outside the interpolation set.

Returns:

L – Two-dimensional NumPy array of shape (N - K, K), where N and K are the numbers of points in the full sampling grid and the Cartesian interpolation grid, respectively. Columns follow the supplied interpolation-index order. Rows correspond to all sampling-grid points outside the Cartesian interpolation grid. Both grids are flattened in NumPy C order, with the last variable varying fastest.

pymor.algorithms.loewner.loewner_quadruple(left_nodes, right_nodes, left_values, right_values, *, left_directions=None, right_directions=None, derivatives=None, force_real=False)

Construct a Loewner quadruple from partitioned transfer function samples.

Assemble the Loewner matrix, shifted Loewner matrix and left and right interpolation data as in [ALI17]. Supports SISO data, full-block MIMO data and tangential MIMO data. No sampling, partitioning or conjugate completion is performed.

Note

In the tangential case, before realification,

\[V_i = \ell_i^T H(\mu_i), \qquad W_j = H(\lambda_j)r_j, \qquad \mathbb{L}_{ij} = \frac{V_i r_j - \ell_i^T W_j}{\mu_i - \lambda_j}.\]

Directions are neither normalised nor conjugated internally. For full-block MIMO data, rows are ordered by left node then output, and columns by right node then input. For a square quadruple, the descriptor-system sign convention is \(E=-L\), \(A=-L_s\), \(B=V\) and \(C=W\).

Parameters:

• left_nodes – Nonempty one-dimensional NumPy array of left interpolation nodes \(\mu_i\), of shape (n_left,).

• right_nodes – Nonempty one-dimensional NumPy array of right interpolation nodes \(\lambda_j\), of shape (n_right,).

• left_values – NumPy array containing transfer function values at left_nodes. Shape (n_left,) for SISO data or (n_left, p, m) for a system with p outputs and m inputs. With tangential directions, already projected values V of shape (n_left, m) are also accepted.

• right_values – NumPy array containing transfer function values at right_nodes. Shape (n_right,) for SISO data or (n_right, p, m) for MIMO data. With tangential directions, already projected values W of shape (p, n_right) are also accepted.

• left_directions – Optional NumPy array of left tangential directions of shape (n_left, p). Requires right_directions. Each row stores a direction \(\ell_i^T\), applied without complex conjugation.

• right_directions – Optional NumPy array of right tangential directions of shape (n_right, m). Requires left_directions. Each row stores a direction \(r_j^T\), used as a column when multiplying a transfer function value.

• derivatives – Optional NumPy array of unprojected transfer function derivatives with respect to the complex argument at left_nodes, with the shape (n_left, p, m) for MIMO data. Required when a left and right node coincide. Entries at other nodes are ignored. Tangential projections of these derivatives are performed internally.

• force_real – If True, transform the quadruple to real arrays using unitary conjugate-pair transformations. Each node set must already contain unique complex conjugate pairs, and values, used derivatives and directions at conjugate nodes must be conjugates. Values and directions at real nodes must be real. No conjugate data are appended.

Returns:

• L – Loewner matrix as a NumPy array. Shape (n_left * p, n_right * m) for full-block MIMO data, or (n_left, n_right) for SISO or tangential data.

• Ls – Shifted Loewner matrix as a NumPy array of the same shape as L.

• V – Left interpolation data as a NumPy array. Shape (n_left * p, m) for full-block MIMO data, (n_left, m) for tangential data, or (n_left, 1) for SISO data.

• W – Right interpolation data as a NumPy array. Shape (p, n_right * m) for full-block MIMO data, (p, n_right) for tangential data, or (1, n_right) for SISO data.

Raises:

AccuracyError – If force_real=True cannot produce real matrices up to roundoff.

pymor.algorithms.loewner.partition_frequencies(nodes, samples, partitioning='even-odd', ordering='regular', force_real=True)

Partition frequency samples into left and right Loewner data sets.

Parameters:

• nodes – One-dimensional NumPy array of complex sampling nodes of shape (n,).

• samples – NumPy array of sampled values of shape (n, ...), aligned with nodes. Used to determine sample norms when ordering is 'magnitude'.

• partitioning –

  Partitioning rule or a tuple (left_indices, right_indices) of index arrays. Available rules are:

  • 'even-odd': assign alternating entries to the left and right sets.

  • 'half-half': assign the first half to the left set and the remainder to the right set. For an odd number of entries, the left set receives one extra entry.

  Explicit index arrays are returned without validation or reordering; ordering and force_real are ignored in this case.

• ordering –

  Ordering applied before a partitioning rule:

  • 'regular': preserve the input order; no sorting by frequency is performed.

  • 'magnitude': order by increasing sample norm.

  • 'random': use a reproducible random ordering with seed zero.

• force_real – If True, keep complex conjugate nodes in the same set. The nodes and samples must already be closed under conjugation. Real nodes and nodes in the upper half-plane are ordered and partitioned separately; the corresponding lower half-plane nodes are then appended to each set. If False, all nodes are ordered and partitioned together.

Returns:

• left_indices – One-dimensional NumPy array of indices into nodes and samples for the left set.

• right_indices – One-dimensional NumPy array of indices into nodes and samples for the right set.

pymor.algorithms.loewner.sample_transfer_function(sampling_values, fom, *, derivative=False)

Sample a TransferFunction or its derivative on a Cartesian grid.

Parameters:

• sampling_values – A one-dimensional NumPy array or a sequence of such arrays. The first array contains Laplace-variable values; subsequent arrays contain parameter values.

• fom – A TransferFunction or a model with a transfer_function attribute.

• derivative – If True, sample the derivative with respect to the complex frequency argument.

Returns:

samples – Sample data of shape tuple(map(len, sampling_values)) + (dim_output, dim_input). A single frequency array produces shape (n, dim_output, dim_input).

pymor.algorithms.loewner.logger