I have a model that has its initial bounding box calculated by finding the maximum and minimum on the x, y and z axes. Producing a correct result like so:

The vertices are then stored in a VBO and only altered with matrices for rotation and bone animation. Currently the bounds are not updated when the model is altered. So the animated and rotated model has bounds like so:

(Maybe it's hard to tell, but the bounds are the same as before, and don't accurately represent the rotated/animated model)

So my question is, how can I calculate the bounding box using the armature matrices and rotation/translation matrices for each model? Keep in mind the modified vertex data is not available because those calculations are performed on the GPU in the shader.

The end result I want is to have an accurate AABB the represents the animated model for picking/basic collision checks.