Efficient way to build a certain square matrix

Question

For odd n I'm looking for a short and swift way to construct with (f.e.)

n = 3

n = 11

n=1111

Answer 1

Here's my take using NestList

cm[n_] := NestList[# + 1 &, Join[Range[n/2 + 1], Reverse@Range[n/2]], n - 1]

Then

cm[11]

Here's a FoldList version (just as fast):

cmf[n_] := FoldList[#1 + #2 &, Join[Range[n/2 + 1], Reverse@Range[n/2]], 
            ConstantArray[1, n - 1]]

The above methods according to the benchmarks posted are already as fast as the fastest methods. Here's a method that's just as fast as Belisarius's (I took a page from his rule based solution):

a4[n_] := With[{t = Join[Range[n/2 + 1], Reverse@Range[n/2]]}, Array[s + # &, n, 0] /. s -> t]

Timings:

I have taken the fastest methods from the five answers with the fastest times and for each function I averaged 7 runs per input for 9 different sizes (1001 to 9001), here are the results:

Edit

After working with JacobAkerbbom to get his LibraryLink method working on Windows see here, I've now included his method which obviously is now the king of the hill.

Answer 2

Edit: See end of post for latest performance enhancement.

f=With[{c = Ceiling[#/2]}, c - 1 + Array[#1 - Abs[c - #2] &, {#, #}]] &;

f[5]

(* {{1, 2, 3, 2, 1}, {2, 3, 4, 3, 2}, {3, 4, 5, 4, 3}, {4, 5, 6, 5, 4}, {5, 6, 7, 6, 5}} *)

Short, sweet, fast.

For more speed,

f5 = With[{c = Ceiling[#/2]},
   Subtract[
    ArrayPad[ConstantArray[Range[#, # + c - 1], c], {{c - 1, 0}, {0, c - 1}}, "Reflected"], 
    Range[# - 1, 0, -1]]] &

And even faster...

f6[n_] := ConstantArray[Join[Range[n/2 + 1], Reverse@Range[n/2]], n] + Range[0, n - 1]

Update: For large n, this seems to do quite well (takes a hit for smallish n before auto-compile kicks in, could just compile overall):

frx[n_] := With[{rng = Range@(n + Floor[n/2]),fl = Floor[n/2], r = Range@n},
   Table[Join[rng[[x ;; x + fl]], rng[[Subtract[x + fl, 1] ;; x ;; -1]]], {x, r}]];

Finally, if this is something you're going to call repeatedly with varying n, particulary larger n, consider injecting a precomputed base of maximum needed size into a lookup function (see later fastest only benchmark for performance illustration):

x2 = With[{base = fr[5001]}, Drop[base, # + 1 ;;, {Ceiling[#/2] + 1, -Ceiling[#/2]}] &];

As requested, a quick bench chart. Usual loungebook caveats apply, and I did not include nasser's or aky's, both get slow on large problems. Looks like f6 and Runnykine`s solutions trade blows for fastest. Update: I'd spaced on Szabolc's answer, added it - seems it and F6 vie for fastest (I'd write off differences to noise), both having a slight but consistent edge over Runnykine's solution (and that itself is small).

Update 2: Added belisarius' second method. Quite quick!

Here's a new set of benchmarks of just the fastest few. Note I've used Timing vs AbsoluteTiming - I think this a much better indicator of efficiency for this kind of test. That said, the earlier ones seemed pretty consistent with Szabolcs' benchmarks, so beyond platform / version / AbsoluteTiming clouding things, hypothesis non fingo why Runnykine's benchmarks seem to differ significantly. Times are average of ten runs per size per tested.

Code I used for this Q-N-D:

tested = {x2,f6, frx, szabolcs, belisaurius2, runnykinea4};
times = 10;
cur = {}
Monitor[
 fastestBmark = 
  Table[cur = {size, (Mean[
         Table[(ClearSystemCache[]; 
           First@Timing@(#@size)), {times}]]) & /@ tested}, {size, 
    301, 5101, 200}], cur]

Finally (really...), one can extend the latter idea if the need is repeated calls over a wide range of n by caching "touchstone" base cases and using the appropriate case for a given n (since dropping has overhead inverse to amount dropped - one could also switch between Drop and Part depending on n vs touchstone size):

x3 = With[{base = 
     frx[#] & /@ {101, 301, 501, 701, 1001, 2001, 3001, 4001, 5001}}, 
   With[{p = Position[{101, 301, 501, 701, 1001, 2001, 3001, 4001, 5001}, 
                       x_ /; x >= #, 1, 1] &@#}, 
     If[p == {}, {}, 
          Drop[base[[p[[1, 1]]]], # + 1 ;;, {Ceiling[#/2] + 1, -Ceiling[#/2]}]]] &];

Here's a quick test against the fastest answer excluding my others:

Answer 3

This is pretty straightforward and very easy to follow even for someone who just started learning Mathematica. This has its value when you need to read your code a year later, even if you're an experienced user.

n = 11;
k = (n + 1)/2;
row = k - Abs[k - Range[n]];
Table[row + i, {i, 0, n - 1}]

Should be fast enough for most application.

---

Benchmarks

The timings from the plot below are outdated. Several people have posted faster versions since, but I haven't had time to update the plot. I'll keep the benchmark code below in case someone wants to experiment with it.

nasser[n_] := Module[{mid},
  mid = Ceiling[n/2];
  SparseArray[{{i_, j_} /; j > mid, {i_, j_} /; j <= mid} :> {n - j + 
      i, i + j - 1}, {n, n}]
  ]

kuba1[n_] :=     
 ArrayPad[HankelMatrix[Range@n, 
    Range[n, n + Floor[n/2]]][[;; , ;; Ceiling[n/2]]], {{0, 0}, {0, 
    Floor[n/2]}}, "Reflected"]

kuba2[n_] := 
 ArrayPad[Array[Range[#, # + Floor[n/2]] &, n], {{0, 0}, {0, Floor[n/2]}}, "Reflected"]


kuba3[n_] := 
 Transpose[Range[n] + # & /@ Join[#, #[[-2 ;; 1 ;; -1]]] &@Range[0, Floor[n/2]]]

rasher1 = 
 With[{c = Ceiling[#/2]}, c - 1 + Array[#1 - Abs[c - #2] &, {#, #}]] &


rasher5 = 
  With[{c = Ceiling[#/2]}, 
    Subtract[
     ArrayPad[
      ConstantArray[Range[#, # + c - 1], c], {{c - 1, 0}, {0, c - 1}},
       "Reflected"], Range[# - 1, 0, -1]]] &;

rasher6[n_] := 
 ConstantArray[Join[Range[n/2 + 1], Reverse@Range[n/2]], n] + 
  Range[0, n - 1]

belisarius[n_] := (Join[#, Rest@Reverse@#] &@Range[n/2 + 1]) + # & /@ Range[0, n - 1]

belisarius2[n_] := With[{k = (Join[#, Rest@Reverse@#] &@Range[n/2 + 1])}, s + Range[0, n - 1] /. s -> k]

algohi[n_] := 
 Module[{k},
  k = Table[i, {i, 1 + #, n + #}] & /@ Range[0, n/2]; 
  Transpose@Join[Most@k, Reverse@k]
 ]

runnykine[n_] := 
 NestList[# + 1 &, Join[Range[n/2 + 1], Reverse@Range[n/2]], n - 1]

szabolcs[n_] := 
 Module[{k, row},
  k = (n + 1)/2;
  row = k - Abs[k - Range[n]];
  Table[row + i, {i, 0, n - 1}]
 ]

Verify they all do what they should:

functions = {"nasser", "kuba1", "kuba2", "kuba3", "rasher1", "rasher5", "rasher6", "belisarius", "belisarius2", "algohi", "runnykine", "szabolcs"};

Equal @@ Through[(ToExpression /@ functions)[11]]
(* True *)

Benchmark, be careful to choose odd n values only:

bench[fun_String] := 
 Module[{f = ToExpression[fun]}, 
  Table[{n, Min@Table[First@AbsoluteTiming[f[n]], {3}]}, {n, 2 Round[2^Range[4, 11, 1/2]] + 1}]
]

results = bench /@ functions;

ListLogLogPlot[results, Joined -> True, 
 PlotRange -> {10^-5, 1}, PlotLegends -> functions, 
 PlotMarkers -> Automatic, PlotStyle -> colours]

Observe how both my and Algohi's solution gets a sudden boost at a certain size threshold. This is due to Table automatically compiling its argument above SystemOptions["CompileOptions" -> "TableCompileLength"].

Further speedups are certainly possible by manually compiling some of the other solutions.

---

Numerical results from my machine, so you don't have to re-run it if you just want to re-plot it.

Caveat: This was run on a laptop. Laptop CPUs have the habit of throttling their frequency when they heat up, which may affect results.

{{{33, 0.006806}, {47, 0.008816}, {65, 0.013106}, {91, 
   0.022375}, {129, 0.044451}, {183, 0.096939}, {257, 0.196829}, {363,
    0.413977}, {513, 0.885362}, {725, 1.889206}, {1025, 
   4.141960}, {1449, 9.229286}, {2049, 20.484070}, {2897, 
   47.915846}, {4097, 123.148882}}, {{33, 0.000029}, {47, 
   0.000046}, {65, 0.000079}, {91, 0.000140}, {129, 0.000270}, {183, 
   0.000531}, {257, 0.001022}, {363, 0.002025}, {513, 0.004224}, {725,
    0.008757}, {1025, 0.019102}, {1449, 0.039411}, {2049, 
   0.078217}, {2897, 0.153589}, {4097, 0.338157}}, {{33, 
   0.000081}, {47, 0.000115}, {65, 0.000159}, {91, 0.000229}, {129, 
   0.000351}, {183, 0.000594}, {257, 0.002159}, {363, 0.004096}, {513,
    0.008323}, {725, 0.016285}, {1025, 0.035317}, {1449, 
   0.070283}, {2049, 0.149849}, {2897, 0.288016}, {4097, 
   0.607098}}, {{33, 0.000096}, {47, 0.000136}, {65, 0.000194}, {91, 
   0.000294}, {129, 0.000477}, {183, 0.000937}, {257, 0.001827}, {363,
    0.003632}, {513, 0.007367}, {725, 0.014935}, {1025, 
   0.034238}, {1449, 0.066544}, {2049, 0.135146}, {2897, 
   0.267405}, {4097, 0.665833}}, {{33, 0.000211}, {47, 0.000402}, {65,
    0.000751}, {91, 0.001488}, {129, 0.002925}, {183, 0.005780}, {257,
    0.011469}, {363, 0.022900}, {513, 0.045602}, {725, 
   0.091659}, {1025, 0.185086}, {1449, 0.370420}, {2049, 
   0.745876}, {2897, 1.493142}, {4097, 3.090481}}, {{33, 
   0.000041}, {47, 0.000070}, {65, 0.000123}, {91, 0.000231}, {129, 
   0.000453}, {183, 0.000896}, {257, 0.001742}, {363, 0.003533}, {513,
    0.006994}, {725, 0.013959}, {1025, 0.029870}, {1449, 
   0.059805}, {2049, 0.120957}, {2897, 0.242442}, {4097, 
   0.567620}}, {{33, 0.000023}, {47, 0.000033}, {65, 0.000054}, {91, 
   0.000090}, {129, 0.000155}, {183, 0.000297}, {257, 0.000567}, {363,
    0.001158}, {513, 0.002287}, {725, 0.005820}, {1025, 
   0.013108}, {1449, 0.025564}, {2049, 0.051499}, {2897, 
   0.100180}, {4097, 0.319559}}, {{33, 0.000198}, {47, 0.000282}, {65,
    0.000398}, {91, 0.000582}, {129, 0.000879}, {183, 0.001349}, {257,
    0.002032}, {363, 0.003323}, {513, 0.005294}, {725, 
   0.008768}, {1025, 0.016870}, {1449, 0.030577}, {2049, 
   0.082289}, {2897, 0.154232}, {4097, 0.285149}}, {{33, 
   0.000086}, {47, 0.000120}, {65, 0.000168}, {91, 0.000249}, {129, 
   0.000392}, {183, 0.000618}, {257, 0.001007}, {363, 0.001706}, {513,
    0.002993}, {725, 0.005540}, {1025, 0.011265}, {1449, 
   0.021607}, {2049, 0.068486}, {2897, 0.128332}, {4097, 
   0.242391}}, {{33, 0.000112}, {47, 0.000175}, {65, 0.000286}, {91, 
   0.000493}, {129, 0.000909}, {183, 0.001677}, {257, 0.001292}, {363,
    0.002553}, {513, 0.005325}, {725, 0.011771}, {1025, 
   0.028106}, {1449, 0.055791}, {2049, 0.110606}, {2897, 
   0.217962}, {4097, 0.534833}}, {{33, 0.000074}, {47, 0.000108}, {65,
    0.000138}, {91, 0.000206}, {129, 0.000188}, {183, 0.000375}, {257,
    0.000625}, {363, 0.001197}, {513, 0.002524}, {725, 
   0.006462}, {1025, 0.014223}, {1449, 0.027980}, {2049, 
   0.056060}, {2897, 0.111271}, {4097, 0.342107}}, {{33, 
   0.000087}, {47, 0.000109}, {65, 0.000152}, {91, 0.000204}, {129, 
   0.000330}, {183, 0.000533}, {257, 0.000695}, {363, 0.001227}, {513,
    0.002388}, {725, 0.006340}, {1025, 0.013557}, {1449, 
   0.025133}, {2049, 0.047835}, {2897, 0.094955}, {4097, 0.262013}}}

Answer 4

n = 11;
mid = Ceiling[n/2];
mat = SparseArray[{{i_, j_} /; j > mid, {i_, j_} /; j <= mid}:>{n-j+i,i+j-1}, {n, n}];

MatrixForm@mat

Answer 5

b[n_] := (Join[#, Rest@Reverse@#] &@Range[n/2 + 1]) + # & /@ Range[0, n-1]
b[11] // MatrixForm

Edit

Enhanced for some speedup, and curiously enough, it competes well with the fastest answers so far (see @rasher's benchamrk):

bs[n_] := With[{k = (Join[#, Rest@Reverse@#] &@Range[n/2 + 1])},  s+Range[0, n-1] /.s-> k]

Answer 6

n = 11;
k = Table[i, {i, 1 + #, n + #}] & /@ Range[0, n/2];
(Transpose@Join[Most@k, Reverse@k]) // TableForm

or

n = 11;
Table[j, {i, 1, 
   n}, {j, (n - 1)/2 - Abs[Range[-(n - 1)/2, (n - 1)/2]] + 
    i}] // TableForm

Answer 7

n = 11;
ArrayPad[
  HankelMatrix[Range@n, Range[n, n + Floor[n/2]]][[;; , ;; Ceiling[n/2]]],
  {{0, 0}, {0, Floor[n/2]}}, 
  "Reflected"]

or

ArrayPad[
  Array[Range[#, # + Floor[n/2]] &, n],
  {{0, 0}, {0, Floor[n/2]}}, 
  "Reflected"]

or

Transpose[
 Range[n] + # & /@ Join[#, #[[-2 ;; 1 ;; -1]]] &@Range[0, Floor[n/2]] 
 ]

Answer 8

Nothing special here, but as there wasn't any solution using Outer, I thought I'd post this:

With[{n = 11}, (* adjust n *)
Outer[#1 + #2 - 1 &, Range[n], Range[n/2 + 1]~Join~Reverse@Range[n/2]]]

Answer 9

Since no answer used the "Extrapolated" option for ArrayPad, which in my opinion is very straightforward (although not as efficient as the top ones), here is my trying:

Clear[extrpol]
extrpol[n_?OddQ] := 
    With[{k = (n + 1)/2}, 
         ArrayPad[2 k - {# + 1, #} & @ {2, 1, 2}, k - {{2, 1}, 2}, "Extrapolated"]
        ]

Answer 10

g[n_] := With[{m=(n+1)/2},Table[n+j-m-Abs[m-i],{i,n},{j,n}]]
g[11]//Transpose//MatrixForm

Answer 11

<< SymbolicC`
<< Developer`
<< CCompilerDriver`
<< CCodeGenerator`

Please don't mind these unnecessary abstractions.

type = "mint";

abstractFunctionName = "makeMatr";

mainFunctionName = abstractFunctionName <> "I_T";

argumentSingletonGetterFunctionName[type_String] := 
  StringJoin["MArgument_get", type];

getter = argumentSingletonGetterFunctionName["Integer"];

typeSpecWL = "MType_Integer";

dataGetterAbstractor[type_String] :=
 "MTensor_get" <> type <> "Data"

dataGetter = dataGetterAbstractor["Integer"];

Generate some SymbolicC

makeMatrSC =
 CFunction[
  "int",
  mainFunctionName,
  {{"WolframLibraryData", "libData"}, {"mint", 
    "Argc"}, {CPointerType["MArgument"], "Args"}, {"MArgument", "Res"}}
  ,
  CBlock[
   {
    CDeclare["int", CAssign["err", "LIBRARY_NO_ERROR"]],
    CDeclare[type, "input"],
    CDeclare["MTensor", "result"], 
    CDeclare[type, CArray["resultDimensions", 2]],

    CAssign[
     "input", 
     CCall[getter, CArray["Args", 0]]
     ],

    CAssign[CArray["resultDimensions", 0], "input"], 
    CAssign[CArray["resultDimensions", 1], "input"], 
    CAssign["err", 
     CCall[CPointerMember["libData", "MTensor_new"], {typeSpecWL, 2, 
       "resultDimensions", CAddress["result"]}]],

    CDeclare[CPointerType[type], "resultDataPtr"],
    CDeclare[type, CAssign["value", 1]],
    CDeclare[type, 
     CAssign["square", COperator[Times, {"input", "input"}]]],
    CDeclare[type, CAssign["half", COperator[Divide, {"input", 2}]]],

    CAssign["resultDataPtr", 
     CCall[CPointerMember["libData", dataGetter], {"result"}]], 
    CDeclare[type, CAssign["iter", 1]],
    CDeclare[type, CAssign["iterRest", 1]],
    CDeclare[type, 
     CAssign["inputMOne" , COperator[Minus, {"input", 1}]]],

    CWhile[
     COperator[LessEqual, {"iter", "square"}],
     CBlock[
      {
       CAssign[CDereference["resultDataPtr"], "value"],
       COperator[Increment, "resultDataPtr"],
       CIf[
        COperator[LessEqual, {"iterRest", "half"}],
        COperator[Increment, "value"],
        COperator[Decrement, "value"]
        ],
       COperator[Increment, "iter"],
       CIf[
        COperator[Equal, {"iterRest", "inputMOne" }],
        CAssign["iterRest", 0],
        COperator[Increment, "iterRest"]
        ]
       }
      ]
     ],
    CCall["MArgument_setMTensor", {"Res", "result"}],
    CReturn["err"]
    }
   ]
  ];

Make it into a string

cCodeString = "DLLEXPORT"<> " " <> ToCCodeString[makeMatrSC];

boilerPlate = "
#include \"WolframLibrary.h\"

/* Return the version of Library Link */
DLLEXPORT mint WolframLibrary_getVersion( ) {
\treturn WolframLibraryVersion;
}

/* Initialize Library */
DLLEXPORT int WolframLibrary_initialize( WolframLibraryData \
libData) {
\treturn LIBRARY_NO_ERROR;
}

/* Uninitialize Library */
DLLEXPORT void WolframLibrary_uninitialize( WolframLibraryData \
libData) {
\treturn;
}

";

totalCString = boilerPlate <> cCodeString;

Create a library (code can be reused, as it makes different versions of the library. Todo: make code that cleans up libraries)

If[! ValueQ[counter], counter = 1;];
counter++;
counterString = ToString[counter];
libraryName = abstractFunctionName <> "Lib" <> counterString;
lib = CreateLibrary[totalCString, libraryName];

Load the library

LibraryLoad[libraryName]

makeMatrLL = 
 LibraryFunctionLoad[libraryName, 
  mainFunctionName, {{Integer}}, {Integer, 2}]

and profit

makeMatrLL[5]

Timing

My previous timings were kind of bad. This is because rI = (rI+1)%input is much slower than

if(rI == input -1){rI = 0}else{rI++}`

So I guess I learned something :).

My code seems to be the fastest code in this Q&A. Below I compare with RunnyKine, but really you should see his answer for a nice graph.

Timing data

timeTable[func_] := 
 Table[Mean@Table[First@Timing@func[in], {5}], {in, 1001, 9001, 1000}]

timings = timeTable[makeMatrLL];
a4Timings = timeTable[a4];
{timings, a4Timings}

 {{0.002772, 0.009899, 0.021801, 0.038221, 0.109890, 0.159009, 
   0.213708, 0.284732, 0.377101}, {0.005047, 0.028429, 0.057490, 
   0.101514, 0.152389, 0.212622, 0.293324, 0.378144, 0.418646}}

Again, see RunnyKine's answer for a nice graph :)

Further improvements

It seems the main thing that is unnecessary in this code is that we check the value of restIter at every iteration. Maybe we can assume that the input size will always be larger than a certain number. In this case, we can inline some code that many times in order to avoid "branches", i.e. using CIf too much. Using a bunch of Ranges in Mathematica also has to deal with conditions when copying the Ranges into the big matrix. My LibraryLink code does not have to copy data, which is something normal Mathematica code cannot avoid. But it seems these conditionals take much more time than copying of data.

Finally Oleksandr made a good point that copying ranges from a single range array would probably be faster.