  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 about 4 million color-selected main sequence turn-off stars. 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 40kpc of (3.7+/-1.2)x10^8 M_sun. The density
profile of the stellar halo is approximately r^{-3}; it is possible that the
power law slope is shallower inside 20kpc and steeper outside that radius. 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 >~100pc,
after accounting for the (known) contribution of Poisson uncertainties. The
fractional RMS deviation 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 analysing the simulations in the same way as the 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 towards a picture in which an important fraction of the
Milky Way's stellar halo has been accreted from satellite galaxies.
