  We define U(n) gauge theory on fuzzy S^2_N x S^2_N as a multi-matrix model,
which reduces to ordinary Yang-Mills theory on S^2 x S^2 in the commutative
limit N -> infinity. The model can be used as a regularization of gauge theory
on noncommutative R^4_\theta in a particular scaling limit, which is studied in
detail. We also find topologically non-trivial U(1) solutions, which reduce to
the known "fluxon" solutions in the limit of R^4_\theta, reproducing their full
moduli space. Other solutions which can be interpreted as 2-dimensional branes
are also found. The quantization of the model is defined non-perturbatively in
terms of a path integral which is finite. A gauge-fixed BRST-invariant action
is given as well. Fermions in the fundamental representation of the gauge group
are included using a formulation based on SO(6), by defining a fuzzy Dirac
operator which reduces to the standard Dirac operator on S^2 x S^2 in the
commutative limit. The chirality operator and Weyl spinors are also introduced.
