# 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)

cost about 0.7 seconds.

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. Commented Jun 28, 2018 at 3:09
• Does using HadamardMatrix[4096, WorkingPrecision -> MachinePrecision] alter the timing significantly? (I’m on mobile and can’t test it myself). Commented Jun 28, 2018 at 3:19
• Also, Mathematica produces two flavors: Method->"BitComplement" and Method->"Sequency". The former is much faster. Which does Matlab produce?
– JimB
Commented Jun 28, 2018 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. Commented Jun 28, 2018 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. Commented Jun 28, 2018 at 5:09

## 2 Answers

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 Commented Jun 28, 2018 at 10:13
• This is a little bit faster: Nest[KroneckerProduct[{{1, 1}, {1, -1}}, #] &, {{1}}, BitLength[4096] - 1] Commented Sep 27, 2018 at 5:30

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.