3
$\begingroup$

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} *)
$\endgroup$
1
  • $\begingroup$ For searching purposes: "skew diagonals" are also called "antidiagonals". $\endgroup$ Sep 21, 2020 at 8:49

3 Answers 3

6
$\begingroup$

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
$\endgroup$
1
  • $\begingroup$ great answer, thank you! $\endgroup$
    – I.M.
    Sep 21, 2020 at 12:24
3
$\begingroup$

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]

enter image description here

enter image description here

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}];

enter image description here

enter image description here

$\endgroup$
1
  • $\begingroup$ great, thank you, I've accepted LouisB's answer since it appeared first. $\endgroup$
    – I.M.
    Sep 21, 2020 at 12:25
0
$\begingroup$

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]}]

enter image description here

$\endgroup$

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge that you have read and understand our privacy policy and code of conduct.

Not the answer you're looking for? Browse other questions tagged or ask your own question.