sliding window on a vector, and multiplying matrix (python)

Question

I would like to multiply sections of a 1D array by a matrix. lets assume a large 1D array size 1xN, and a nxn matrix (such as n<N). If I take adjacent sections size 1xn, the result is again 1xN array (id N\n is an integer) however, I would like to perform a sliding window operation on this array. Again, the operation is matrix multiplication. which means that in every step of the operation I multiply a section of the vector (size 1xn) by nxn matrix. the question is how to summarize the results to get again 1xN array. I assume that since every cell of the original array is being used n-1 times, I need to average all of these results.

Graphical representation, with n=3: enter image description here Now we need to perform row-averaging.

enter image description here the question is - how do I build such a 2D array for the averaging?

---------------- suggestion---------------------

One can look at it as some kind of convolution, where instead of a dot product, we use a matrix product. In this "convolution", the input is a 1D array, the kernel is a (square) matrix, and the operation is matrix multiplication.

More specifically:

every step of the operation takes a 1Xn section of the input array, and multiplies it with the kernel matrix, to provide a 1Xn result.

the sliding operation is trivial. the summation is a question.

Answer 1

import numpy as np
from scipy.signal import oaconvolve
import matplotlib.pyplot as plt

# Create a random 1xN series
in1 = np.random.randn(1000).cumsum()

# Create two filter kernels (n x n, with `n`=3) for demonstration:
# 1. Copy kernel
kernel_copy = np.array([[0, 0, 0], [0, 1, 0], [0, 0, 0]])
# Example: Median over 3 points
kernel_med3 = np.array([[1, 0, 0], [0, 1, 0], [0, 0, 1]])

# Nomalizing, so that the sum of each kernel is 1
kernel_copy = kernel_copy / kernel_copy.sum()
kernel_med3 = kernel_med3 / kernel_med3.sum()

# Pad time series to avoid discontinuities at the edges
# (Filter assumes zeros, which leads to unwanted "step function")
# Edge values will be repeated for 3 points
pad = kernel_copy.shape[0]  # `n`=3
in1_padded = np.pad(in1, pad, "edge")

# Doing overlap-and-add convolution (zero padded FFT convolve)
# `mode` must be "full", so that point-wise sum over axis 1 works correctly.
# `oaconvolve` calculates n x 1D-convolutions. The resulting shape is n x (N+2*padding-1) which must finally point-wise aggregated by summation to 1 x (N+2*padding-1)
# The result is finally sliced, to remove the previous padded edges
y1 = oaconvolve(in1_padded[:, None], kernel_copy, mode="full").sum(axis=1)[pad+1:-(pad+1)]
y2 = oaconvolve(in1_padded[:, None], kernel_med3, mode="full").sum(axis=1)[pad+1:-(pad+1)]

# Plot results
plt.plot(in1, 'k-', label="original")
plt.plot(y1, 'g--', label="copy")
plt.plot(y2, 'b--', label="med3")
plt.legend()
plt.show()

Answer 2

If you don't care whether result is a 3×4 full matrix, or a 6×4 matrix with holes (as your answer to my comment seems to suggest), then, what you describe is just a matrix multiplication with every extract from the vector.

So what you need is an array whose columns are all extracts (1st column is x₀, x₁, x₂; 2nd column is x₁, x₂, x₃; ...). And then you just have to multiply you matrix by this one, to get a 3×4 matrix whose columns are M×[x₀,x₁,x₂], M×[x₁,x₂,x₃]...

To create such an array, there is sliding_window_view function in np.lib.stride_tricks,

Let start with a minimal reproducible example

N=3
M=np.arange(N*N).reshape(-1,N) # Just a matrix
# array([[0, 1, 2],
#        [3, 4, 5],
#        [6, 7, 8]])
X=np.array([5,6,7,8,9,10]) # Just a vector

Then

np.lib.stride_tricks.sliding_window_view(X,3)

is

array([[ 5.,  6.,  7.],
       [ 6.,  7.,  8.],
       [ 7.,  8.,  9.],
       [ 8.,  9., 10.]])

Which is the array I wanted (extract of size 3 window from X), but in rows rather than columns. It is easy from there to transpose it and multiply M by it.

Note that this is not really an array, but a view to another one. So this occupy no memory in RAM. It is just another view on X. Data are those of X. Memory usage is the one of X.

M @ np.lib.stride_tricks.sliding_window_view(X,3).T

is

array([[ 20.,  23.,  26.,  29.],
       [ 74.,  86.,  98., 110.],
       [128., 149., 170., 191.]])

So that is my one-liner

One-liner to produce a 3×4 matrix

M @ np.lib.stride_tricks.sliding_window_view(X,N).T

To produce a 6×4 matrix

If you really want a 6x4 matrix (your answer seems to indicate that you don't care, and even that a 3×4 matrix is better, since your objective seems to be averaging that — but you seemed unsure, yet did not say much about your real objective. On another hand, all you said is that 6×4 is OK, not whether 3×4 also is or not), then we have to rearrange that result.

Note that if we look at your drawing, and compare to that result, and reading them from left to right and then from top to bottom (so like Chinese writing if I am not mistaken), the difference between your drawing and my result is that in yours you have 4 (6-3+1) between the last element of a column and the first of the next.

So lets add those 4 0

np.pad(M @ np.lib.stride_tricks.sliding_window_view(X,3).T, ((0,4),(0,0)), constant_values=np.nan)

is

array([[ 20.,  23.,  26.,  29.],
       [ 74.,  86.,  98., 110.],
       [128., 149., 170., 191.],
       [ nan,  nan,  nan,  nan],
       [ nan,  nan,  nan,  nan],
       [ nan,  nan,  nan,  nan],
       [ nan,  nan,  nan,  nan]])

Now if we could read this left to right, top to bottom, and retranscribe it in a 6×4 matrix, we would have what we want. We can't do that. But resize does it for rows (in latin reading, from top to bottom and left to right). So all we have to do is to transpose that, and resize it

np.resize(np.pad(M @ np.lib.stride_tricks.sliding_window_view(X,3).T, ((0,4),(0,0)), constant_values=np.nan).T, (4,6))

array([[ 20.,  74., 128.,  nan,  nan,  nan],
       [ nan,  23.,  86., 149.,  nan,  nan],
       [ nan,  nan,  26.,  98., 170.,  nan],
       [ nan,  nan,  nan,  29., 110., 191.]])

Almost there. Since we have transposed before resize, we need to transpose after also

np.resize(np.pad(M @ np.lib.stride_tricks.sliding_window_view(X,3).T, ((0,4),(0,0)), constant_values=np.nan).T, (4,6)).T

array([[ 20.,  nan,  nan,  nan],
       [ 74.,  23.,  nan,  nan],
       [128.,  86.,  26.,  nan],
       [ nan, 149.,  98.,  29.],
       [ nan,  nan, 170., 110.],
       [ nan,  nan,  nan, 191.]])

Or, with generic dimensions

np.resize(np.pad(M @ np.lib.stride_tricks.sliding_window_view(X,3).T, ((0,len(X)-N+1),(0,0)), constant_values=np.nan).T, (len(X)-N+1,len(X))).T

tl;dr

My two one-liners are

M @ np.lib.stride_tricks.sliding_window_view(X,N).T

If you want a full matrix, with no empty cells

np.resize(np.pad(M @ np.lib.stride_tricks.sliding_window_view(X,3).T, ((0,len(X)-N+1),(0,0)), constant_values=np.nan).T, (len(X)-N+1,len(X))).T

If you really want that version with sliding result columns and useless nan values.

Results are the same, but for the nan padding.

If you have a choice, 1st one is simpler, faster, use less memory.