Abstract 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 Tr ( A ρ ( t ) ) in the ensemble thermalize: For N → ∞ every such function tends to the value Tr ( A ρ eq ( ∞ ) ) + Tr ( A ρ ( 0 ) ) g 2 ( t ) . Here ρ ( t ) is the time-dependent density matrix of the system, A is a Hermitean operator standing for an observable, and ρ eq ( ∞ ) 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 | 3 / 2 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 of random matrices. Comparison with the ‘eigenstate thermalization hypothesis’ of (Srednicki 1999 J. Phys. A: Math. Gen. 32 1163) shows overall agreement but raises significant questions.