CODEWITHSUNDEEP
algorithmmatrixoptimizationprocessing-efficiencytoeplitz
Harry
Harry

Reputation: 938

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

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

enter image description here

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}")

Upvotes: 1

Views: 44

Answers (1)

Reinderien
Reinderien

Reputation: 15231

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))

Upvotes: 0

Related Questions