Loewner Reductor

[lbalicki]

This PR implements a reductor based on the Loewner interpolation framework. This is a draft for now because there might be some room for discussion regarding:

1. _partition_frequencies is more or less an attempt to generalize the different partitioning choices considered here. In the chapter it seems like frequencies are considered to be complex only (makes sense, although it would be good to at least take 0 into consideration) and also they don't seem to care about keeping complex conjugate pairs together which is something I tried to do in the implementation.
2. There is the loewner_quadruple function for all related matrices. However, it maybe makes sense to have separate functions to generate the individual matrices. As I couldn't think of a use case I left them together for now.
3. There is also room for optimization in the setup of the matrices.

[artpelling]

Great! What are your thoughts regarding derivative data?

[lbalicki]

Great! What are your thoughts regarding derivative data?

Do you mean as in TFBHIReductor? It is not really part of the "basic" Loewner, right? My thought was to get the Loewner reductor into pyMOR and then do a refactoring PR organizing Loewner related things from TFBHIReductor, PAAAReductor and LoewnerReductor. But of course I would be open to approach this differently.

[codecov]

Codecov Report

Merging #1952 (fef2c19) into main (f42f279) will decrease coverage by 0.07%.
The diff coverage is 77.89%.

❗ Current head fef2c19 differs from pull request most recent head 18fc887. Consider uploading reports for the commit 18fc887 to get more accurate results

Additional details and impacted files

Flag	Coverage Δ
github_actions	76.14% <77.89%> (+0.01%)	⬆️
gitlab_ci	87.53% <77.89%> (-0.07%)	⬇️

Flags with carried forward coverage won't be shown. Click here to find out more.

Impacted Files	Coverage Δ
src/pymortests/demos.py	92.48% <ø> (ø)
src/pymor/reductors/loewner.py	74.21% <74.21%> (ø)
src/pymordemos/dd_heat.py	96.77% <96.77%> (ø)

... and 1 file with indirect coverage changes

[artpelling]

Nice job so far. I've added some comments.

[artpelling]

The SVDs should be cached. Especially in the current implementation where the partitioning and data does not change. It might also be nice to pass the partitioning parameter to the reduce method, but then a more sophisticated caching mechanism would be needed.

[lbalicki]

I added caching just now. I will have to think about passing the partitioning to reduce. It makes sense, but would also mean having some caching mechanism that maps partitionings to corresponding SVDs... Maybe it's easier to leave it as is.

[artpelling]

Nice work. Thanks for adding the full MIMO version! I like the mimo_handling variable :)

[artpelling]

@lbalicki I think there still is an error in the construction of the Loewner quadruple for mimo_handling='full'. V seems to be correct, but L, Ls, W have the correct entries, but seem to be permuted somehow.
In particular, the blocks of W are not equal to the transfer function values.

This is the code snippet that I am using which is correct as far as I can tell:

    if shifted:
        L = (np.expand_dims(s[i], axis=(1, 2)) * Hs[i])[:, np.newaxis] - (np.expand_dims(s[j], axis=(1, 2)) * Hs[j])[np.newaxis]
    else:
        L = Hs[i][:, np.newaxis] - Hs[j][np.newaxis]

    L /= np.expand_dims(np.subtract.outer(s[i], s[j]), axis=(2, 3))

    return np.concatenate(np.concatenate(L, axis=1), axis=1)

Also, the spectra of the Loewner matrices are different...

[artpelling]

Sorry, I've misclicked. Did not mean to close the PR.

Here is a MWE that you can check the old vs. new behaviour and see the difference.

#!/usr/bin/env python3

from pymor.models.transfer_function import TransferFunction
import numpy as np
import matplotlib.pyplot as plt
from pymor.models.iosys import LTIModel
from pymor.reductors.loewner import LoewnerReductor, _partition_frequencies


def _construct_loewner_matrix(s, Hs, dHs=None, shifted=False, partitioning='even-odd', ordering='regular'):
    """Constructs a (shifted) Loewner matrix as a |NumPy array|."""
    assert isinstance(s, np.ndarray)
    if hasattr(Hs, 'transfer_function'):
        Hs = Hs.transfer_function
    assert isinstance(Hs, TransferFunction) or isinstance(Hs, np.ndarray) and Hs.shape[0] == s.shape[0]
    assert dHs is None or isinstance(dHs, np.ndarray)

    assert isinstance(shifted, bool)
    assert partitioning in ('even-odd', 'half-half', 'same') or len(partitioning) == 2
    if partitioning == 'same':
        raise NotImplementedError
    if isinstance(partitioning, str):
        partitioning = _partition_frequencies(s, Hs, partitioning=partitioning, ordering=ordering)
    assert ordering in ('magnitude', 'random', 'regular')

    # compute derivatives if needed
    if not isinstance(partitioning, str) and np.any(np.isin(*partitioning)) or partitioning == 'same':
        if isinstance(Hs, TransferFunction):
            assert Hs.dtf is not None
            dHs = np.stack([Hs.eval_dtf(si) for si in s])
        else:
            assert isinstance(dHs, np.ndarray) and Hs.shape == dHs.shape

    if isinstance(Hs, TransferFunction):
        Hs = np.stack([Hs.eval_tf(si) for si in s])

    i, j = _partition_frequencies(s, Hs, partitioning, ordering) if isinstance(partitioning, str) else partitioning

    if shifted:
        L = (np.expand_dims(s[i], axis=(1, 2)) * Hs[i])[:, np.newaxis] - (np.expand_dims(s[j], axis=(1, 2)) * Hs[j])[np.newaxis]
    else:
        L = Hs[i][:, np.newaxis] - Hs[j][np.newaxis]

    L /= np.expand_dims(np.subtract.outer(s[i], s[j]), axis=(2, 3))

    return np.concatenate(np.concatenate(L, axis=1), axis=1)


def loewner_quadruple(s, sys, partitioning='even-odd', ordering='regular'):
    part = _partition_frequencies(s, sys, partitioning=partitioning, ordering=ordering)
    IL = _construct_loewner_matrix(s, sys, partitioning=part)
    ILs = _construct_loewner_matrix(s, sys, partitioning=part, shifted=True)
    IV = np.vstack([sys.eval_tf(si) for si in s[part[0]]])
    IW = np.hstack([sys.eval_tf(si) for si in s[part[1]]])
    return IL, ILs, IV, IW


r = 10
n = 100
s = np.geomspace(1, 100, n)*1j
fom = LTIModel.from_matrices(np.random.rand(n, n), np.random.rand(n,3), np.random.rand(4,n))
lr = LoewnerReductor(s, fom, mimo_handling='full', conjugate=False)
rom = lr.reduce(r=r)

L, Ls, V, W = loewner_quadruple(s, fom.transfer_function)
LhLs = np.hstack([L, Ls])
Y, S1, _ = spla.svd(LhLs, full_matrices=False)
LvLs = np.vstack([L, Ls])
_, S2, Xh = spla.svd(LvLs, full_matrices=False)
Yhr = Y[:, :r].conj().T
Xr = Xh[:r, :].conj().T
rom2 = LTIModel.from_matrices(- Yhr @ Ls @ Xr, Yhr @ V, W @ Xr, D=None, E=- Yhr @ L @ Xr)


wlim = (s[0]/1j, s[-1]/1j)
fig1, ax = plt.subplots()
fom.transfer_function.mag_plot(wlim, ax=ax)
rom.transfer_function.mag_plot(wlim, ax=ax, label='reductor')
rom2.transfer_function.mag_plot(wlim, ax=ax, label='old-code')
plt.legend()

fig2, ax = plt.subplots()
ax.semilogy(lr.loewner_svds[1], c='orange', label='reductor')
ax.semilogy(lr.loewner_svds[2], c='orange', label='reductor')
ax.semilogy(S1, c='g', label='old-code')
ax.semilogy(S2, c='g', label='old-code')
plt.legend()
plt.show()

[pmli]

Below are some minor comments.

[lbalicki]

Thanks for the comments. @artpelling thanks for providing the code for the MIMO case, it seems like there was an issue with my reshaping in the end. Also, the conversion to real Loewner matrices caused some issues. I think it works now, but I still need think about how to make the reshaping code look nicer.

[pmli]

I think it would be good to test the different options. Running

pytest --cov --cov-config=setup.cfg --cov-report=html src/pymortests/demos.py
open htmlcov/index.html

on my machine says that 57% of the code is covered.

[lbalicki]

I added a simple test for the LoewnerReductor similar to the one for MTReductor. Not sure if this will increase the coverage in src/pymortests/demos.py but most of the code should be covered in the tests with this.

[pmli]

I added a simple test for the LoewnerReductor similar to the one for MTReductor. Not sure if this will increase the coverage in src/pymortests/demos.py but most of the code should be covered in the tests with this.

Thanks. The coverage is 90% when running

pytest --cov --cov-config=setup.cfg --cov-report=html src/pymortests/reductors/loewner.py