How does one practically calculate the non-linearity of a multi-output boolean function like the AES s-box? I am working on multi-output boolean function i.e. a function that takes an $n$-bit input and gives an $m$-bit output for all possible inputs (i.e. $2^n$ inputs).

To be more exact: I am trying to do a cryptanalysis of the AES 8x8 s-box, with the goal to learn the whole procedure of finding/calculating non-linearity. 

Doing some research, I learned that I need to find the minimum hamming distance from all $n$ variable affine functions. Sadly enough, I don't  really grasp all of the math which would be required to find the non-linearity. 

In the end, I’ld like to be able to calculate the Linear Approximation Table of the AES s-box. Can someone please explain how to build a LAT of the AES Sbox. How can I practically calculate the non-linearity of a multi-output boolean function like the AES s-box?
Before we start with vectorial Boolean functions, let's recall the definition of the nonlinearity of a Boolean function:

$$\mathcal{NL}(f) = \min_{a \in \mathbb{F}_2^n} d_H(f, \ell_a \oplus b),$$

where $\ell_a \oplus b$ represents the affine Boolean function defined by the bitvector $a$: $\ell_a(x) = a \cdot x$ ($\cdot$ is the dot product).
The above equation pretty much defines the nonlinearity of a Boolean function as the minimum Hamming distance $d_H$ to some affine function.

This distance $d_H$ can be expressed using the Walsh-Spectrum of $f$.
That is, $d_H(f, \ell_a) = 2^{n - 1} - \frac{1}{2}\mathcal{W}_f(a)$.
To see why this is true, the definition of the Walsh transform should help:
$$\mathcal{W}_f(a) = \sum_{x \in \mathbb{F}_2^n} (-1)^{f(x) \oplus \ell_a(x)}.$$
(and note that $(-1)^{g(x)} = 1 - 2g(x)$, $\sum_x g(x) = w_H(x))$)

Minimizing the distance hence corresponds to maximizing the Walsh-Spectrum:
$$\mathcal{NL}(f) = 2^{n - 1} - \frac{1}{2} \max_{a \in \mathbb{F}_2^n} \left|\mathcal{W_f(a)}\right|.$$

Why the absolute value? Remember that we're measuring the distance to affine functions and note that $d_H(f, \ell_a \oplus 1) = 2^{n - 1} + \frac{1}{2}\mathcal{W}_f(a)$.

Why is this useful? Well, it turns out there is a rather efficient algorithm to compute the Walsh-Spectrum of some Boolean function. See for example here.

Now that we can compute the nonlinearity of a Boolean function, it's easy to define that of a vectorial Boolean function $F$:

$$\mathcal{NL}(F) = \min_{a \in \mathbb{F}_2^n} \mathcal{NL}(F \cdot a).$$

In other words, the nonlinearity is the minimum of the nonlinearities of the Boolean functions which are linear combinations of the coordinates (outputs) of $F$.

Given what we know about the nonlinearity of Boolean functions, we can compute this as

$$\mathcal{NL}(F) =  2^{n - 1} - \frac{1}{2} \max_{a \in \mathbb{F}_2^n, b \in \mathbb{F}_2^m} \left|\mathcal{W}_{F\cdot b}(a)\right|.$$

Sometimes the expression in the absolute value is called the Walsh transform of $F$:
$$\mathcal{W}_F(a, b) = \mathcal{W}_{F\cdot b}(a)$$

To compute this, you can use the fast Walsh-Hadamard transform (FWHT) algorithm again.
Note also this corresponds directly to the LAT (depending on conventions, the term $2^{n - 1}$ is added or not).