# Mean values of skew diagonals of a $(n+1,n)$ matrix

How mean values of skew diagonals of a $$(n+1,n)$$ matrix can be computed efficiently?

Here is my naive implementation:

ClearAll[build] ;
build[matrix_] := Block[
{col,row,signal},
{col,row} = Dimensions[matrix] ;
signal = ConstantArray[0,2*row] ;
Do[
signal[[i]] = Table[If[q+p==i+1,matrix[[q,p]],Nothing],{q,1,col},{p,1,row}] ;
signal[[i]] = Mean[Flatten[signal[[i]]]] ;
,{i,1,2*row,1}
] ;
signal
]


Looks like it's time complexity is $$O(n^3)$$, can it be reduced?

Example:

n = 4 ;
ncols = n + 1 ;
nrows = n ;
matrix = Array[m,{ncols,nrows}] ;
matrix
build[matrix]
(* {{m[1,1],m[1,2],m[1,3],m[1,4]},{m[2,1],m[2,2],m[2,3],m[2,4]},{m[3,1],m[3,2],m[3,3],m[3,4]},{m[4,1],m[4,2],m[4,3],m[4,4]},{m[5,1],m[5,2],m[5,3],m[5,4]}} *)
(* {m[1,1],1/2 (m[1,2]+m[2,1]),1/3 (m[1,3]+m[2,2]+m[3,1]),1/4 (m[1,4]+m[2,3]+m[3,2]+m[4,1]),1/4 (m[2,4]+m[3,3]+m[4,2]+m[5,1]),1/3 (m[3,4]+m[4,3]+m[5,2]),1/2 (m[4,4]+m[5,3]),m[5,4]} *)

n = 4 ;
ncols = n + 1 ;
nrows = n ;
data = Range[1,2*n] ;
data = Partition[data,n,1] ;
data
build[data]
(* {{1,2,3,4},{2,3,4,5},{3,4,5,6},{4,5,6,7},{5,6,7,8}} *)
(* {1,2,3,4,5,6,7,8} *)

• For searching purposes: "skew diagonals" are also called "antidiagonals". Commented Sep 21, 2020 at 8:49

We can use the Diagonal function, but first we must "rotate" the matrix. Start by constructing a matrix with $$n$$ rows and $$n+1$$ columns:

n = 3;
mat = Table[m[irow, jcol], {irow, n}, {jcol, n + 1}];
mat // MatrixForm


$$\left( \begin{array}{cccc} m(1,1) & m(1,2) & m(1,3) & m(1,4) \\ m(2,1) & m(2,2) & m(2,3) & m(2,4) \\ m(3,1) & m(3,2) & m(3,3) & m(3,4) \\ \end{array} \right)$$

"Rotate" the elements like this

rot = Reverse @ Transpose @ mat;
rot // MatrixForm


$$\left( \begin{array}{ccc} m(1,4) & m(2,4) & m(3,4) \\ m(1,3) & m(2,3) & m(3,3) \\ m(1,2) & m(2,2) & m(3,2) \\ m(1,1) & m(2,1) & m(3,1) \\ \end{array} \right)$$

The diagonals of rot are the skew diagonals of mat. So the mean of the skew diagonals of mat can be obtained by

diags = Table[Diagonal[rot, k], {k, 1-n,n-2];
Mean /@ diags  // Column


$$\begin{array}{l} \frac{1}{2} (m(1,2)+m(2,1)) \\ \frac{1}{3} (m(1,3)+m(2,2)+m(3,1)) \\ \frac{1}{3} (m(1,4)+m(2,3)+m(3,2)) \\ \frac{1}{2} (m(2,4)+m(3,3)) \\ \end{array}$$

For $$n+1$$ rows and $$n$$ columns in the original matrix, use

n = 3;
mat = Table[m[irow, jcol], {irow, n + 1}, {jcol, n}];
rot = Reverse @ Transpose @ mat;
diags = Table[Diagonal[rot, k], {k, 2 - n, n - 1}];
Mean /@ diags

• great answer, thank you!
– I.M.
Commented Sep 21, 2020 at 12:24

The following is long to comment.

Numeric: power 2 behavior

t0 = Table[H = RandomReal[{-1, 1}, {n, n + 1}];
{n, AbsoluteTiming[HR = Reverse /@ H;
Table[
Total@Diagonal[HR, i], {i, -Length@HR +
1, +Length@HR}];][[1]]}, {n, 1000, 10000, 1000}];
ff = a x^n /. FindFit[t0, a x^n, {a, n}, x]


array: around power 2 behavior

t0 = Table[H = Array[aa, {n, n + 1}];
{n, AbsoluteTiming[HR = Reverse /@ H;
Table[
Total@Diagonal[HR, i], {i, -Length@HR +
1, +Length@HR}];][[1]]}, {n, 100, 1000, 100}];


• great, thank you, I've accepted LouisB's answer since it appeared first.
– I.M.
Commented Sep 21, 2020 at 12:25

Here is a procedural style solution.

ClearAll[build] ;
build[matrix_] := Block[
{row,col,array, signal, shift, start, count, i, j},
{row,col} = Dimensions[matrix] ;
array = Flatten[Transpose[matrix]]    ;
signal = Table[0,2*col]  ;
shift = 1 ;
Do[
{
start = i ;
If[
i > row,
start = i + col*shift ;
shift++ ;
] ;
count = 0 ;
Do[
{
signal[[i]] += array[[start+j*col]] ;
count++ ;
},
{j,0,Min[i,col]-shift,1}
] ;
signal[[i]]/=count ;
},
{i,1,2*col,1}
] ;
signal
] ;


Compare scaling with Diagonal solution:

ClearAll[LouisB] ;
LouisB[matrix_] := Block[
{row, col, rotate, count},
{row, col} = Dimensions[matrix] ;
rotate = Reverse[Transpose[matrix]] ;
Table[Mean[Diagonal[rotate,count]],{count,1-col,row-1}]
] ;

ClearAll[test] ;
test[function_,seed_:"SkewSum"] := Block[
{data, fit},
SeedRandom[seed] ;
data = Table[
Block[
{matrix, time},
matrix = RandomReal[{-1, 1}, {n+1, n}] ;
time = First[AbsoluteTiming[function[matrix]]] ;
{n, time}
],
{n, 500, 5000, 500}
] ;
fit = NonlinearModelFit[data,a*x^n,{a,n},x] ;
fit = fit["Function"] ;
Show[
ListPlot[data, PlotTheme -> "Detailed", ImageSize -> 400, PlotStyle -> Directive[{PointSize[Large], Red}], AspectRatio -> 1/2],
Plot[fit[x],{x,500,10000}, PlotStyle -> Black],
PlotLabel -> fit[x]
]
] ;

Row[{test[LouisB],test[build]}]