Think about what the inertia tensor is doing- it describes an objects resistance to rotational motion along each of the 3 axes. Each of the 3 elements which comprise the diagonal of the inertia tensor describe the resistance to motion along each of the 3 local space axes. If you made some of the elements in the column vectors non-zero, you would simply be limiting one of the other two axes- which is the natural result of applying a rotation to the inertia tensor.
Since any given object must have an origin and original orientation, it only makes sense to define them in local space as simply as possible, for both memory reasons as well as simplicity.
I suggest reviewing the polyhedral mass properties algorithm to gain a better understanding of how the diagonal is calculated. It is essentially a volume integral of a solid defined by triangular regions.
I should also mention that the inverse inertia tensor is the quantity actually used by a physics engine, so defining them as a diagonal makes the inverse process significantly simpler.