We study the operator norm discrepancy of i.i.d. random matrices, initiating the matrix-valued analog of a long line of work on the $\ell^{\infty}$ norm discrepancy of i.i.d. random vectors. First, using repurposed results on vector discrepancy and new first moment method calculations, we give upper and lower bounds on the discrepancy of random matrices. We treat i.i.d. matrices drawn from the Gaussian orthogonal ensemble (GOE) and low-rank Gaussian Wishart distributions. In both cases, for what turns out to be the "critical" number of $Θ(n^2)$ matrices of dimension $n \times n$, we identify the discrepancy up to constant factors. Second, we give a new analysis of the matrix hyperbolic cosine algorithm of Zouzias (2011), a matrix version of an online vector discrepancy algorithm of Spencer (1977) studied for average-case inputs by Bansal and Spencer (2020), for the case of i.i.d. random matrix inputs. We both give a general analysis and extract concrete bounds on the discrepancy achieved by this algorithm for matrices with independent entries (including GOE matrices) and Gaussian Wishart matrices.