  An isolated quantum system is said to thermalize if ${\rm Tr} (A \rho(t)) \to
{\rm Tr} (A \rho_{\rm eq})$ for time $t \to \infty$. Here $\rho(t)$ is the
time-dependent density matrix of the system, $\rho_{\rm eq}$ is the
time-independent density matrix that describes statistical equilibrium, and $A$
is a Hermitean operator standing for an observable. We show that for a system
governed by a random-matrix Hamiltonian (a member of the time-reversal
invariant Gaussian Orthogonal Ensemble (GOE) of random matrices of dimension
$N$), all functions ${\rm Tr} (A \rho(t))$ in the ensemble thermalize: For $N
\to \infty$ every such function tends to the value ${\rm Tr} (A \rho_{\rm
eq}(\infty)) + {\rm Tr} (A \rho(0)) g^2(t)$. Here $\rho_{\rm
  eq}(\infty)$ is the equilibrium density matrix at infinite temperature. The
oscillatory function $g(t)$ is the Fourier transform of the average GOE level
density and falls off as $1 / |t|$ for large $t$. With $g(t) = g(-t)$,
thermalization is symmetric in time. Analogous results, including the symmetry
in time of thermalization, are derived for the time-reversal non-invariant
Gaussian Unitary Ensemble (GUE) of random matrices. Comparison with the
``eigenstate thermalization hypothesis'' of Ref.~\cite{Sre99} shows overall
agreement but raises significant questions.
