`amen_mv` exceeds sweep limit when working with complex numbers

[levdik]

amen_mv weirdly exceeds number of sweeps limit when working with some complex-valued matrices and vectors.

Here is the simplest example I could find:

d = 3

x = tt.rand(2, d)
x += 0.1j * x

a = tt.matrix(tt.rand(4, d))
y = amen_mv(a, x, 1e-12, nswp=1000)[0]

error = np.linalg.norm(y.full(asvector=True) - np.dot(a.full(), x.full(asvector=True)))
print(f'Error: {error}')

This results in:

amen-mv: swp=1{1}, max_dx=1.564e+00, max_r=4
amen-mv: swp=1{0}, max_dx=3.960e-01, max_r=2
amen-mv: swp=2{1}, max_dx=2.027e-01, max_r=4
 
 .
 .
 .

amen-mv: swp=999{0}, max_dx=3.928e-01, max_r=2
amen-mv: swp=1000{1}, max_dx=2.005e-01, max_r=4
amen-mv: swp=1000{0}, max_dx=3.921e-01, max_r=2
Error: 5.999547960677571e-14

The method converges quickly, though, as the error is negligible even if nswp=1:

amen-mv: swp=1{1}, max_dx=3.305e+00, max_r=4
amen-mv: swp=1{0}, max_dx=3.960e-01, max_r=2
Error: 1.1873159181672063e-14

[daskol]

First of all, it worth to note that tolerance tol is not an error of matvec approximation. It characterizes an L₂-norm of change in core during iterations.

We can see this if we run amen_mv() with verb=1 to see norm of residuals during iterations. Or we can see change it with your slightly modified code.

import numpy as np
from numpy.testing import assert_array_almost_equal

import tt
from tt.amen import amen_mv


def test_amen_mv(nosweeps, d=3):
    x = tt.rand(2, d)
    x += 0.1j * x

    a = tt.matrix(tt.rand(4, d))
    y = amen_mv(a, x, 1e-12, nswp=nosweeps, verb=0)[0]

    y_exact = np.dot(a.full(), x.full(asvector=True))
    error = np.linalg.norm(y.full(asvector=True) - y_exact)
    print(f'error: {error:e}')
    return y

Then we set manually random seed and compute matvecs for 10 sweeps and 11 sweeps.

np.random.seed(42)
y_curr = test_amen_mv(10)

np.random.seed(42)
y_next = test_amen_mv(11)

We do not see difference in dense output. Resulting vector are the same.

assert_array_almost_equal(y_curr.full(), y_next.full())  # Ok

However, cores could change significant.

def get_cores(self: tt.vector) -> np.ndarray:
    """Tuple of views on TT-cores. Each element in the list is a view on
    underlying buffer.
    """
    offset = 0
    cores = []
    ranks = tuple(self.r.tolist())
    for core_shape in zip(ranks[:-1], self.n, ranks[1:]):
        core_size = np.prod(core_shape)
        core = self.core[offset:offset + core_size]
        cores.append(core.reshape(core_shape, order='F'))
        offset += core_size
    return cores


cores_curr = get_cores(y_curr)[::-1]
cores_next = get_cores(y_next)[::-1]
for i, (core_curr, core_next) in enumerate(zip(cores_curr, cores_next)):
    print(f'core_curr[{i}] =\n')
    print(core_curr.squeeze())
    print()
    print(f'core_next[{i}] =\n')
    print(core_next.squeeze())
    print()
    assert_array_almost_equal(core_curr, core_next)

The absolute difference in this example is about 1.4 and relative is about 1.5, i.e. elements are significantly differ. This is what are the first TT-cores.

core_curr[0] =

[[-2.49107418-2.74883544e-17j  0.92247455+1.80844437e-18j]
 [ 0.33472097-4.18401217e-02j  0.90388919-1.12986149e-01j]]

core_next[0] =

[[-2.49107418e+00+4.96719386e-17j  9.22474549e-01-5.99197900e-17j]
[ 8.43940704e-19+3.37325845e-01j  5.00939089e-17+9.10923456e-01j]]

This is quite usual situation with ALS-like method in complex arithmetics, I guess. We probably can add another stopping criterion.

UPD May be @kharyuk has some thoughts.

[levdik]

@daskol thank you for the quick response! Now it's clear why such behavior occurs.

Another stopping criterion would be a great improvement, in my opinion. Otherwise, the method's behavior isn't quite intuitive. And, in my case, nswp has to be manually changed constantly.