Speed up the generation of large Hadamard matrices

The MATLAB matrix is similar to the BitComplement type. Avoiding HadamardMatrix makes the construction slightly faster.

r1 = Nest[ArrayFlatten[{{#, #}, {#, -#}}] &, {{1}}, Log2[4096]]; // AbsoluteTiming
r2 = Sign[HadamardMatrix[4096, Method -> "BitComplement",
                         WorkingPrecision -> MachinePrecision]]; // AbsoluteTiming
r1 === r2

{0.14940275, Null}
{0.36388055, Null}
True

Regarding the creation of an image, you should not use MatrixPlot. It constructs a Graphics which is awfully slow. Rather you can create an Image

AbsoluteTiming[
 r1 = Nest[ArrayFlatten[{{#, #}, {#, -#}}] &, {{1}}, Log2[4096]];
 Image[r1]
]

Some spelunking of the code for HadamardMatrix[] reveals why it is so slow with the option setting Method -> "Sequency". (Coolwater's answer already shows the method corresponding to Method -> "BitComplement").

Altho the function is internally using a CompiledFunction[] for the purpose, it is slow because of repeated calls to IntegerDigits[]. I have thus reorganized the algorithm being used in the code that follows:

hadmat[n_Integer?Positive, prec_: ∞] /; BitAnd[n, n - 1] == 0 := 
   With[{bits = IntegerDigits[Range[0, n - 1], 2, BitLength[n] - 1]},
        ((-1)^(ListCorrelate[{{1, 1}}, Reverse[bits, 2], 1, {0}].Transpose[bits]))/
        Sqrt[N[n, prec]]]

where I have used the classical bit-twiddling test for a power of 2.

Even if it is uncompiled, hadmat[] handily outperforms the built-in:

AbsoluteTiming[t1 = hadmat[4096];]
   {4.54179, Null}

AbsoluteTiming[t2 = HadamardMatrix[4096];]
   {37.9819, Null}

t1 === t2
   True

(It should perhaps be emphasized that the two method settings for HadamardMatrix[] give different sets of matrices, but both produce matrices with the Hadamard property.)

---

For giggles, I wrote a compiled routine for generating the Hadamard matrix:

hmc = Compile[{{n, _Integer}, {l, _Integer}},
              With[{bits = Table[IntegerDigits[k - 1, 2, l], {k, n}]},
                   (-1)^(bits.Table[Switch[l - j - k + 1, 0, 1, 1, 1, _, 0],
                                           {j, l}, {k, l}].Transpose[bits])]];

hadm2[n_Integer?Positive, prec_: ∞] /; BitAnd[n, n - 1] == 0 := 
    hmc[n, BitLength[n] - 1]/Sqrt[N[n, prec]]

where convolution is replaced with multiplication by an antitriangular matrix. As expected, this method is slower than the method based on ListCorrelate[], but is still faster than the built-in.