How can I write the Anderson hamiltonian as a matrix? How can I write this Hamiltonian:
$$ H = \sum E_d \hat{n}_d + \sum_k \epsilon_k\hat{n}_k + \sum_k V_{kd} (\hat{a}^\dagger_k \hat{a}_d + \hat{a}^\dagger_d \hat{a}_k) $$
in matrix form using its eigenvector representation. $ V_{kd} $ are coupling constants, $ \hat{a}^\dagger $ and $ \hat{a} $ are Fermion creation and annihilation operators, and $ \hat{n} = \hat{a}^\dagger\hat{a} $ is the number operator. 
 A: This Hamiltonian is a second-quantized Hamiltonian. This means that we have a freedom to talk about arbitrary number of particles. If we have this freedom, is obvious that we have a "big" hilbert space (Fock Space). Is not convenient to view this Hamiltonian in a matrix picture because the hamiltonian commute with number operator. If we do that we came with a lot of redundancy. But the first-quantized version of this Hamiltonian is very convenient in the matrix picture.
In this picture we have as the diagonal elements of the matrix the energy of the impurity state and the conduction band states. The off-diagonal elements is zero except for the first line and column. This non-zero off-diagonal elements is responsible to transition of band states with impurity state.
Nevertheless, note that exist a non-bilinear term that you forgot, i.e., a term that can't be written as:
$$
\sum_{i\,j}f_{i}^{\dagger}U_{i}{}_{j}f_{j}
$$
The term is
$$
U\,f_\uparrow ^\dagger f_\uparrow f_\downarrow ^\dagger f_\downarrow
$$
so, this matrix picture can't appreciate the correct Anderson hamiltonian. What happens is that this term is crucial for the Kondo effect as well as all the many-body interesting phenomena.
The second quantization matrix is very complicated. If we use the symmetries of this hamiltonian as charge conservation and spin conservation we can construct a infinity number of manageable matrices, by work on subspace of a fixed charge and spin. So we can construct the real matrix as a block diagonal matrix like this:
$$
\left(\begin{array}{cccc}
0 & 0 & 0 & 0\\
0 & \left(\begin{array}{cc}
\varepsilon_{d} & J\\
J & \textrm{conduction band}
\end{array}\right) & 0 & 0\\
0 & 0 & \left(\textrm{Wigner-Eckart sum of spins}\right) & 0\\
0 & 0 & 0 & \ddots
\end{array}\right)
$$
The first matrix is actually doubly because the two possible spins $\pm \frac{1}{2}$. This matrix represent the one-particle matrix, and is very easy and non exciting. For construct the other matrix in a manageable way we need to use the Wigner-Eckart theorem and be careful with the non-bilinear term (i.s. the coulomb term). In the limit $N\rightarrow \infty$ is very difficult to work on this matrices.
Exist a suitable discretization (logarithm discretization) of the band conduction that yields numerical tractable way to work on in this matrices. Furthermore, work on this discretization provide to us a basis that tri-diagonalize the conduction part in such a way that only one state of this discretization couples direct with the impurity, giving an approximate hamiltonian:
$$
H =H_{\textrm{impurity}}+\overline{D}\left(c_{d}^{\dagger}f_{0}+f_{0}^{\dagger}c_{d}\right)+D\sum_{N=0}^{\infty}\Lambda^{-\frac{N}{2}}\left(f_{N+1}^{\dagger}f_{N}+f_{N}^{\dagger}f_{N+1}\right)
$$
This procedure are called numerical renormalization group, and the key idea is to control the resolution in energy scale of your model by truncating the matrix. Also, we can use this procedure to maping the renormalization flow. And finally, work iteratively to diagonalize the matrices produced by the possible truncations.
