It is probably most intuitive to rotate the object around the axis perpendicular to the current drag direction, either incrementally with each mouse motion, or relative to the drag start position. The two options give slightly different user interactions, which each have their pluses and minuses.
There is a relatively straightforward way to convert an angle and a 3d vector representing the axis being rotated around into a rotation matrix.
You are right in that updating a raw rotation matrix through incremental rotations will result in the matrix no longer being a pure rotation matrix. That is because a 3x3 rotation matrix has three times as much data than needed to represent a rotation.
A more compact way to represent rotations is with Euler Angles, having a minimal 3 value vector. You could take the current rotation as an Euler angle vector, convert it to a matrix, apply the rotation (incremental or otherwise), and convert the matrix back to an Euler angle vector. That last step would naturally eliminate any non-rotational component to your matrix, so that you once again end up with a pure rotational matrix for the next state.
Euler angles are conceptually nice, however it is a lot of work to do the back and forth conversions.
A more practical choice is Quaternions (also), which are four element vectors. The four elements specify rotation and uniform scale, and it happens that if you go in and normalize the vector to unit length, you will get a scale factor of 1.0. It turns out that an angle-axis value can also be converted to a quaternion value very easily by
q.x = sin(0.5*angle) * axis.x;
q.y = sin(0.5*angle) * axis.y;
q.z = sin(0.5*angle) * axis.z;
q.w = cos(0.5*angle);
You can then take the quaternion product (which uses only simple multiplication and addition) of the current rotation quaternion and incremental rotation quaternion to get a new quaternion which represents performing both rotations. At that point you can normalize the length to ensure a pure rotation, but otherwise continue iteratively combining rotations.
Converting the quaternion to a rotation matrix is very straightforward (uses only multiplication and addition) when you want to display the model in its rotated state using traditional graphics API's.