We give a systematic study of the Hardy spaces of functions with values in the non-commutative Lp-spaces associated with a semifinite von Neumann algebra M. This is motivated by matrix valued harmonic analysis (operator weighted norm inequalities, operator Hilbert transform), as well as by the recent development of non-commutative martingale inequalities. Our non-commutative Hardy spaces are defined by non-commutative Lusin integral functions. It is proved in this dissertation that they are equivalent to those defined by the non-commutative Littlewood-Paley G-functions. We also study the Lp boundedness of operator valued dyadic paraproducts and prove that their Lq boundedness implies their Lp boundedness for all 1 < q < p < ∞.