Efficient way to sum the columns of a Toeplitz/Hankel matrix

Question

There are many repeated values in an (N × N) Toeplitz matrix. Example:

The same is true for its mirror image, the Hankel matrix.

I am trying to find the most efficient way to sum the columns of such matrices, where I define "most efficient" to mean the minimal number of add/subtract operations.

I have included my Python attempt for the Hankel matrix below. If I have counted correctly, this needs 3N-4 additions. Is it possible to do better?

Clarification: This Python code is only meant as (functioning) pseudocode. I'm not targeting a CPU, so the execution speed of this code doesn't matter - only the number of add/sub operations.

import numpy as np import scipy # NxN matrix dimension N = 8 # Define the last row and first column of the NxN Hankel matrix r = np.random.randint(1, 100, N) # Last row c = np.concatenate((np.random.randint(1, 100, N-1), [r[0]])) # First column # Cumulative sums cumsum_r = np.cumsum(r[1:]) cumsum_c = np.cumsum(np.flip(c)) # Compute column sums sums = np.zeros(N) for i in range(N): if i == 0: sums[i] = cumsum_c[N-1] else: sums[i] = cumsum_c[N-1-i] + cumsum_r[i-1] # Explicitly construct the Hankel matrix and check the sums A = scipy.linalg.hankel(c, r) mismatch = np.sum(A, axis=0) - sums print(f"Mismatch: {mismatch}")

The lower-bound for the number of addition needed is at least 2N-1 since each different value needs to be read once. — Jérôme Richard
– Jérôme Richard, Commented Nov 13, 2024 at 22:25
You can do that by just reading the first line and column: in 2N-1 additions plus some multiplications (2N-3 multiplications by constants). a*5+b*4+c*3+... so something like np.sum(line * np.arange(line.size, 0, -1)) + np.sum(column[1:] * np.arange(column.size-1, 0, -1)) — Jérôme Richard
– Jérôme Richard, Commented Nov 13, 2024 at 22:29
An alternative solution is to compute partial sums so to avoid some subtraction and get something between 2N and 3N add/sub without any multiplications. Maybe a Fenwick tree can help to do that. — Jérôme Richard
– Jérôme Richard, Commented Nov 13, 2024 at 22:33

Reinderien · Accepted Answer · 2024-11-13 02:14:34Z

Seems reasonable, but don't loop, and assert your mismatch rather than printing it:

import numpy as np import scipy # NxN matrix dimension N = 8 # Define the last row and first column of the NxN Hankel matrix rand = np.random.default_rng(seed=0) r = rand.integers(low=1, high=100, size=N) # Last row c = np.concatenate(( # First column rand.integers(low=1, high=100, size=N-1), r[:1], )) # Cumulative sums cumsum_r = np.cumsum(r[1:]) cumsum_c = np.cumsum(c[::-1]) # Compute column sums sums = np.empty(N) sums[0] = cumsum_c[N-1] sums[1:] = cumsum_c[-2::-1] + cumsum_r # Explicitly construct the Hankel matrix and check the sums A = scipy.linalg.hankel(c, r) assert np.array_equal(sums, np.sum(A, axis=0))

Thanks for the tips regarding my horrible Python coding. If I understand correctly, this still uses the same number of additions. I added a clarification to my question to note that my Python code was only intended as pseudocode to illustrate my attempt.

Collectives™ on Stack Overflow

Efficient way to sum the columns of a Toeplitz/Hankel matrix

1 Answer 1

1 Comment

Hot Network Questions

Collectives™ on Stack Overflow

1 Answer 1

1 Comment

Related