In preparation for the era of the time-domain astronomy with upcoming large-scale surveys, we propose a state-space representation of a multivariate damped random walk process as a tool to analyze irregularly-spaced multi-filter light curves with heteroscedastic measurement errors. We adopt a computationally efficient and scalable Kalman-filtering approach to evaluate the likelihood function, leading to maximum O(k3n) complexity, where k is the number of available bands and n is the number of unique observation times across the k bands. This is a significant computational advantage over a commonly used univariate Gaussian process that can stack up all multi-band light curves in one vector with maximum O(k3n3) complexity. Using such efficient likelihood computation, we provide both maximum likelihood estimates and Bayesian posterior samples of the model parameters. Three numerical illustrations are presented in this work: (i) analyzing simulated five-band light curves for a comparison with independent single-band fits; (ii) analyzing five-band light curves of a quasar obtained from the Sloan Digital Sky Survey (SDSS) Stripe 82 to estimate short-term variability and timescale; (iii) analyzing gravitationally lensed g- and r-band light curves of Q0957+561 to infer the time delay. Two R packages, Rdrw and timedelay, are publicly available to fit the proposed models.