Implement the 2D Hadamard Transform

Question

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.

Answer 1

Haskell, 111 102 97 92 bytes

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

Answer 2

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\$.

Answer 3

MATL, 12 bytes

&nKYLw/Gy,Y*

Try it online! Or verify all test cases.

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)

Answer 4

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!

Answer 5

APL (Dyalog), 28 bytes

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

Try it online!

Answer 6

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 add each pair  (adds the 2 subvectors)
             -ᵐ              #         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!

Answer 7

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
                $  - last two links as a monad:
               L   -   length (number of rows)
              ÷    -   divide (vectorises)
                 ⁺ - repeat the previous link (divide by the length again)

Answer 8

R, 80 78 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.

Answer 9

Pari/GP, 51 bytes

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

Try it online!