# Implement the 2D Hadamard Transform

The Hadamard transform works similarly to the Fourier transform: it takes a vector and maps it to its frequency components, which are the Walsh functions. Instead of sine waves in the Fourier transform, the Walsh functions are discrete "square waves" so to speak. This makes it easier to compute.

The size of a Hadamard matrix is a power of two, 2^n x 2^n. We can build a Hadamard matrix for a given n recursively. $$H_1 = \begin{pmatrix}1&1\\1 & -1\end{pmatrix}$$ And then build block matrices in the same pattern to get $$H_n = \begin{pmatrix}H_{n-1}&H_{n-1}\\H_{n-1} & -H_{n-1}\end{pmatrix}$$

For instance,

$$H_3 = \begin{pmatrix}1&1&1&1&1&1&1&1\\1&-1&1&-1&1&-1&1&-1\\1&1&-1&-1&1&1&-1&-1\\1&-1&-1&1&1&-1&-1&1\\1&1&1&1&-1&-1&-1&-1\\1&-1&1&-1&-1&1&-1&1\\1&1&-1&-1&-1&-1&1&1\\1&-1&-1&1&-1&1&1&-1\end{pmatrix} = \begin{pmatrix}H_2&H_2\\H_2&-H_2\end{pmatrix}$$

Note how the entries are always 1 or -1 and symmetric! Easy! Then the Hadamard transform is simply multiplying H_n by a column vector and dividing by 2^n.

In order to apply the Hadamard transform to 2D images, we have to apply the 1D HT across every column and then every row. In short, for an input matrix X, the output will be a new matrix Y given by two matrix multiplications:

$$Y = \frac{1}{2^{2n}}H_nXH_n$$

Another way to implement this is by first applying the HT to every column in the image separately. Then apply it to every row, pretending they were column vectors. For example, to compute the 2D HT of $$\begin{pmatrix}2 & 3\\ 2 & 5\end{pmatrix}$$

we first do

$$\frac{1}{2}\begin{pmatrix}1&1\\1 & -1\end{pmatrix}\begin{pmatrix}2 \\ 2\end{pmatrix} = \begin{pmatrix}2 \\ 0\end{pmatrix} \quad\quad\quad \frac{1}{2}\begin{pmatrix}1&1\\1 & -1\end{pmatrix}\begin{pmatrix}3 \\ 5\end{pmatrix} = \begin{pmatrix}4 \\ -1\end{pmatrix} \quad\rightarrow\quad \begin{pmatrix}2&4\\0 & -1\end{pmatrix}$$

which takes care of the columns, and then we do the rows,

$$\frac{1}{2}\begin{pmatrix}1&1\\1 & -1\end{pmatrix}\begin{pmatrix}2 \\ 4\end{pmatrix} = \begin{pmatrix}3 \\ -1\end{pmatrix} \quad\quad\quad \frac{1}{2}\begin{pmatrix}1&1\\1 & -1\end{pmatrix}\begin{pmatrix}0 \\ -1\end{pmatrix} = \begin{pmatrix}-0.5 \\ 0.5\end{pmatrix} \quad\rightarrow\quad \begin{pmatrix}3&-1\\-0.5 & 0.5\end{pmatrix}$$

Finally, one can also compute this recursively using the fast Hadamard transform. It's up to you to find out if this is shorter or not.

## Challenge

Your challenge is to perform the 2D Hadamard transform on a 2D array of arbitrary size 2^n x 2^n, where n > 0. Your code must accept the array as input (2D array or a flattened 1D array). Then perform the 2D HT and output the result.

In the spirit of flexibility, you have a choice for the output. Before multiplying by the prefactor 1/(2^(2n)), the output array is all integers. You may either multiply by the prefactor and output the answer as floats (no formatting necessary), or you may output a string "1/X*" followed by the array of integers, where X is whatever the factor should be. The last test case shows this as an example.

You may not import the Hadamard matrices from some convenient file that already contains them; you must build them inside the program.

## Test Cases

Input:

2,3
2,5


Output:

3,-1
-0.5,0.5


Input:

1,1,1,1
1,1,1,1
1,1,1,1
1,1,1,1


Output:

1,0,0,0
0,0,0,0
0,0,0,0
0,0,0,0


Input:

1,0,0,0
1,1,0,0
1,1,1,0
1,1,1,1


Output:

0.6250,0.1250,0.2500,0
-0.1250,0.1250,0,0
-0.2500,0,0.1250,0.1250
0,0,-0.1250,0.1250


Input:

0,0,0,0,0,0,0,1
0,0,0,0,0,0,2,0
0,0,0,0,0,3,0,0
0,0,0,0,4,0,0,0
0,0,0,5,0,0,0,0
0,0,6,0,0,0,0,0
0,7,0,0,0,0,0,0
8,0,0,0,0,0,0,0


Output:

1/64*
36,4,8,0,16,0,0,0
-4,-36,0,-8,0,-16,0,0
-8,0,-36,-4,0,0,-16,0
0,8,4,36,0,0,0,16
-16,0,0,0,-36,-4,-8,0
0,16,0,0,4,36,0,8
0,0,16,0,8,0,36,4
0,0,0,-16,0,-8,-4,-36


## Scoring

This is code golf, so the smallest program in bytes wins.

• Nice first challenge! – Arnauld Oct 13 '19 at 22:39
• – Luis Mendo Oct 14 '19 at 10:38
• @Arnauld Thanks! I don't think I'm good enough to golf yet, but it's still fun seeing the other answers. – HiddenBabel Oct 14 '19 at 18:13
• "You may not import the Hadamard matrices from some convenient file that already contains them; you must build them inside the program." Does this include using build-in functions that return the matrix? – Kroppeb Oct 14 '19 at 19:12
• an upvote and applause for not requiring floating point – ngn Oct 14 '19 at 19:27

import Data.Matrix
t x|n<-nrows x=h n*x*h n
h 1=1
h i|m<-h$div i 2=(/2)<$>m<->m<|>(m<->(-m))


This uses the Data.Matrix module. The input matrix has to be of type Matrix Float or Matrix Double.

TIO doesn't provide Data.Matrix, so no link.

Edit: -9 bytes thanks to @ceased to turn counterclockwis, -5 bytes thanks to @H.PWiz, -5 bytes thanks to @xnor

• I am amazed that identity 1 is the shortest way to represent $\left[1\right]$, but it does seem to be. – Ad Hoc Garf Hunter Oct 14 '19 at 3:32
• @SriotchilismO'Zaic actually no, pure 1 is shorter. Heck, even 1 does the trick. – ceased to turn counterclockwis Oct 14 '19 at 15:28
• t x|Right m<-inverse$h$nrows x=m*x*m – H.PWiz Oct 14 '19 at 23:45
• Can you divide the matrix by 2 in the recursive step to avoid needing to scale or invert at the end? – xnor Oct 15 '19 at 22:33

# JavaScript (ES6), 114 bytes

This version is golfed a bit further by embedding the function $$\h\$$ within the callback of the deepest loop.

M=>(g=z=>M=M.map((r,y)=>r.map((_,x)=>r.map(v=>s+=(h=m=>m?-h(m&~-m):z?v:M[i-1][x])(i++&(z?x:y)),i=s=0)&&s/i)))(g())


Try it online!

# JavaScript (ES6),  161 ... 117  116 bytes

M=>(h=m=>m?-h(m&~-m):1,g=z=>M=M.map((r,y)=>r.map((_,x)=>r.map(v=>s+=(z?v:M[i][x])*h(i++&(z?x:y)),s=i=0)&&s/i)))(g())


Try it online!

### How?

Each 0-indexed cell $$\(x,y)\$$ of a Hadamard matrix $$\H_n\$$ is equal to:

$$H_n(x,y)=\cases{ 1,&\text{if (x\And y) is evil}\\ -1,&\text{otherwise} }$$

or:

$$H_n(x,y)=(-1)^{\operatorname{popcount}(x \And y)}$$

where $$\\And\$$ is a bitwise AND and $$\\operatorname{popcount(m)}\$$ is the number of bits set in $$\m\$$.

The function $$\h\$$ computes $$\(-1)^{\operatorname{popcount}(m)}\$$ recursively:

h = m => m ? -h(m & ~-m) : 1


(This is similar to my answer to Is this number evil?)

The function $$\g\$$ computes either $$\\frac{1}{2^n}H_nM\$$ or $$\\frac{1}{2^n}MH_n\$$, depending on the flag $$\z\$$:

g = z =>                      // z = flag cleared on the 1st pass, set on the 2nd
M =                         // update M[]:
M.map((r, y) =>           //   for each row r[] at position y in M[]:
r.map((_, x) =>         //     for each cell at position x in r[]:
r.map(v =>            //       for each value v in r[]:
s +=                //         add to s:
(z ? v            //           v = M(i, y) if z is set,
: M[i][x]) *   //             or M(x, i) otherwise
h(i++ & (z ? x    //           multiplied by H(x, i) if z is set,
: y)), //             or H(i, y) otherwise
//           and increment i afterwards
s = i = 0           //         start with s = i = 0
)                     //       end of inner map() over the columns
&& s / i              //       yield s divided by the width of the matrix
)                       //     end of outer map() over the columns
)                         //   end of map() over the rows


We invoke it twice to compute $$\\frac{1}{2^{2n}}H_nMH_n\$$.

• "If (x&y) is evil" :D – Jonathan Allan Oct 14 '19 at 16:31

# MATL, 12 bytes

&nKYLw/Gy,Y*


### Explanation

Consider input [2,3; 2,5] as an example.

&n     % Input (implicit). Push number of rows and number of columns
% STACK: 2, 2
KYL    % Hadamard matrix of specified size
% STACK: 2, [1 1; 1 -1]
w/     % Swap, divide element-wise
% STACK: [0.5 0.5; 0.5 -0.5]
G      % Push input again
% STACK: [0.5 0.5; 0.5 -0.5], [2,3; 2,5]
y      % Duplicate from below
% STACK: [0.5 0.5; 0.5 -0.5], [2,3; 2,5], [0.5 0.5; 0.5 -0.5]
,      % Do twice
Y*   %   Matrix multiply
% End (implicit)
% STACK: [3 -1; -0.5 0.5]
% Display (implicit)

• I don't know MATL. But it might be shorter to get an inverse, instead of dividing by the number of rows. (both get the same result) – H.PWiz Oct 15 '19 at 0:26
• @H.PWiz Thanks for the idea! I didn't know that. Unfortunately, computing the inverse matrix (and getting rid of the extra number of rows) would require more characters – Luis Mendo Oct 15 '19 at 8:47

# Octave, 93 79 bytes

@(x)(r=f(f=@(f)@(x){@()[x=f(f)(x-1) x;x -x],1}{1+~x}())(log2(r=rows(x)))/r)*x*r


This computes $$\\frac{1}{2^n}H_nX\frac{1}{2^n}H_n\$$

Try it online!

# Octave, 47 33 bytes (uses builtin)

@(x)(r=hadamard(r=rows(x))/r)*x*r


Try it online!

# APL (Dyalog), 28 bytes

{(⍵∘⌹+.×⌹)(,⍨⍪⊢,-)⍣(2⍟≢⍵)⍪1}


Try it online!

# Brachylog, 23 bytes

{ḍ\⟨+ᵐ↰₁ᵐ-ᵐ⟩c/₂ᵐ|}ᵐ\↰₁ᵐ

• Input: A list of columns
• Output: A list of rows

## How:

It uses the Fast Hadamard Transformation. AFAIK Brachylog doesn't have matrix multiplication, nor binary operations so I doubt using the "evil numbers technique" will help shorten the code.

    {ḍ\⟨+ᵐ↰₁ᵐ-ᵐ⟩c/₂ᵐ|}ᵐ\↰₁ᵐ   #
{ḍ\⟨+ᵐ↰₁ᵐ-ᵐ⟩c/₂ᵐ|}ᵐ       # Fast transform:
{               }ᵐ       #   for each column:
ḍ                       #     Split in 2
|         #       (if this fails, return the collum.
#         This is the case in the final recursion step)
\                      #      Zip,
-ᵐ              #         and sub each pair  (subs the 2 subvectors)
⟨  ↰₁ᵐ  ⟩              #      and apply the transformation on both results
c             #      concat both transformed results.
/₂^m        #        and half each element
\      #  Transpose
↰₁ᵐ   #    and apply the transformation again



Try it online!

# Jelly, 18 bytes

J’&þB§-*ðæ×⁺@÷L$⁺  A monadic Link accepting a list of lists of integers (any numbers really), the square input matrix of side $$\2^n\$$, which yields a list of lists of floats. Try it online! ### How? It feels like one should be able to outer-product the dot-product of the binary representations of the 0-indexed coordinates, however to do so one would need to left-pad the binary representations with zeros which, I believe, costs too many bytes. J’&þB§-*ðæ×⁺@÷L$⁺ - Link: list of lists of integers, M
J                  - range of length = [1,2,...,rows(M)]
’                 - decrement (vectorises) = [0,1,...,rows(M)-1]
              - use left argument as both arguments of:
þ               -   outer-product with:
&                -     bit-wise AND         [[0&0, 0&1, ...],[1&0, 1&1, ...], ...]
B             - convert to binary (vectorises)
§            - sums (vectorises at depth 1) - i.e. bit-sums
-           - literal minus one
*          - exponentiate (vectorises) - i.e. the appropriately sized Hadamard matrix, H
ð        - start a new dyadic chain - i.e. f(H, M)
æ×       - matrix multiply - i.e. H×M
@     - with swapped arguments i.e. f(H×M, H):
⁺      -   repeat previous link i.e. H×M×H
L   -   length (number of rows)
÷    -   divide (vectorises)
⁺ - repeat the previous link (divide by the length again)


# R, 8078 74 bytes

function(M){for(i in 1:log2(nrow(M)))T=T%x%matrix(1-2*!3:0,2)/2
T%*%M%*%T}


Try it online!

Calculates $$\\frac{1}{2^n}H_nX\frac{1}{2^n}H_n\$$; constructs $$\\frac{1}{2^n}H_n\$$ by Kronecker product %x%.

Looks like I can still golf, recalling that 1-2*!3:0 is shorter than c(1,1,1,-1) as I realized in the comments to this answer.

# Pari/GP, 51 bytes

a->h=1;while(#h<#a,h=matconcat([h,h;h,-h])/2);h*a*h
`

Try it online!