We have used data from the Sloan Digital Sky Survey (SDSS) Data Release 5 to explore the overall structure and substructure of the stellar halo of the Milky Way using ~4 million color-selected main-sequence turnoff stars with 0.2 < g − r < 0.4 and 18.5 ≤ r < 22.5. We fit oblate and triaxial broken power law models to the data, and found a "best-fit" oblateness of the stellar halo 0.5 < c/a < 0.8, and halo stellar masses between galactocentric radii of 1 and 40 kpc of 3.7 ± 1.2 × 108 M☉. The density profile of the stellar halo is approximately ρ r−α, where –2 > α > − 4. Yet, we found that all smooth and symmetric models were very poor fits to the distribution of stellar halo stars because the data exhibit a great deal of spatial substructure. We quantified deviations from a smooth oblate/triaxial model using the rms of the data around the model profile on scales 100 pc, after accounting for the (known) contribution of Poisson uncertainties. Within the DR5 area of the SDSS, the fractional rms deviation σ/total of the actual stellar distribution from any smooth, parameterized halo model is 40%: hence, the stellar halo is highly structured. We compared the observations with simulations of galactic stellar halos formed entirely from the accretion of satellites in a cosmological context by analyzing the simulations in the same way as the SDSS data. While the masses, overall profiles, and degree of substructure in the simulated stellar halos show considerable scatter, the properties and degree of substructure in the Milky Way's halo match well the properties of a "typical" stellar halo built exclusively out of the debris from disrupted satellite galaxies. Our results therefore point toward a picture in which an important fraction of the stellar halo of the Milky Way has been accreted from satellite galaxies.