# Speed up the generation of large Hadamard matrices

I need to generate a Hadamard matrix with dimensions of $$4096$$, but it seems Mathematica is much slower than MATLAB when generating such matrices.

This program

(MatrixPlot @ HadamardMatrix[4096]) // AbsoluteTiming


cost about 36 seconds on my computer to finish. While for MATLAB,

 H = hadamard(4096); imshow(H)


I tested the pure process of generating the H matrix. The code is listed below.

Mathematica cost 20 seconds:

HadamardMatrix[4096, WorkingPrecision -> MachinePrecision] // AbsoluteTiming


MATLAB cost 0.12 second:

hadamard(4096)


I was wondering if there is any method to increase the speed or efficiency of Mathematica on this problem?

Jim has provided a way of changing the Method:

HadamardMatrix[4096, Method -> "BitComplement"] // AbsoluteTiming
HadamardMatrix[4096, Method -> "Sequency"] // AbsoluteTiming


It improves the speed a little bit, but still much slower than MATLAB

• While I don’t think that this is the cause of the slow down, you still shouldn’t include the plotting function in your timing. Secondarily, consider also that the Hadamard matrix produced by Mathematica is normalized by a factor of $1/\sqrt{n}$, whereas Matlab produces a matrix of $\pm 1$. Having said that, if the timing in both systems is accurate, the difference is very large indeed. – MarcoB Jun 28 '18 at 3:09
• Does using HadamardMatrix[4096, WorkingPrecision -> MachinePrecision] alter the timing significantly? (I’m on mobile and can’t test it myself). – MarcoB Jun 28 '18 at 3:19
• Also, Mathematica produces two flavors: Method->"BitComplement" and Method->"Sequency". The former is much faster. Which does Matlab produce? – JimB Jun 28 '18 at 3:23
• @MarcoB，@JimB，It should be right that it should be blamed on the Sqrt[n] here. Method->"BitComplement" and Method->"Sequency" can improve this problem a little bit. However,it seems WorkingPrecision -> MachinePrecision cann't improve the behavoir too much. I' wondering whether there is a solution to force the internal calculation to use machine precision numbers, just as Matlab does. – yulinlinyu Jun 28 '18 at 4:55
• @MarcoB change the working precision is not a very effective way. It still uses an awfully long time to produce the result. – BNHSX Jun 28 '18 at 5:09

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


• That's much faster! Thanks – BNHSX Jun 28 '18 at 10:13
• This is a little bit faster: Nest[KroneckerProduct[{{1, 1}, {1, -1}}, #] &, {{1}}, BitLength[4096] - 1] – J. M. is away Sep 27 '18 at 5:30

Some spelunking of the code for HadamardMatrix[] reveals why it is so slow for the 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 redundant computations to IntegerDigits[]. I have the 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}

{37.9819, Null}

t1 === t2
True


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 is slower than the method based on ListCorrelate[], but this method is still faster than the built-in.