Meaning and Origin of an Expression which Involves Virtual Displacement As an additional point of confusion related to the answer given here:
Confusion with Virtual Displacement 
I have encountered the following expression in my study of virtual displacements.
$$\delta{\vec{r_i}}=\sum_j\frac{\partial\vec{r_i}}{\partial{q_j}}\delta{q_j}.$$
It is difficult for me to even parse this statement. On the left hand side we have an infinitesimal variation at $\vec{r_i}$ (which I now understand to be a tangent vector to the position manifold $M$ at that point) and it is being equated to the sum on the right. The partial derivative in the sum seems easy enough to understand; I assume that it denotes the usual component-wise partial derivative of $\vec{r_i}$ with respect to the coordinate $q_j$. The derivative then seems to be multiplied by an infinitesimal virtual displacement (which I now understand to be a tangent vector to the constraint submanifold $C$). The sum is obviously over all of the coordinates.
My confusion, then, is two-fold. First, what kind of multiplication is being carried out between $\delta{q_j}$ and $\frac{\partial\vec{r_i}}{\partial{q_j}}$? They don't even seem to be in the same spaces. Even if we took the image of $\delta{q_j}$ under some kind of inclusion map, it still isn't clear what operation we are performing. Is it the cross-product, the inner product, or something else?
Lastly, it is unclear to me where this expression even came from. In class, my professor practically pulled it out of a hat with the statement that "it is easy to see". Well, I fail to see how this is in any way obvious. In CuriousMind's comment here:
Mathematics of the Virtual Displacement 
It is mentioned that this is actually the definition of an infinitesimal virtual displacement. However, I fail to see how this "definition" corresponds to my understanding that an infinitesimal virtual displacement  is a tangent vector to a manifold. Any answer or clarification will be much appreciated.
 A: Before you talk about dq or dr you need to talk about q and r.  
These are coordinates on the manifolds.  A typical use of the mechanism you have cited is found in expressing differential elements in curvilinear coordinates.  The {q} are a set of independent parameters and {r} the usual Cartesian coordinates in R^(N).  It should be noted, because I have seen a lot of confusion on this point, that r is NOT a vector, and likewise {q} is not a vector in any space.  These are coordinates and are to be thought of as a set of scalar functions {r_i}, i = 1...N.  When dealing with Euclidean space it is easy to develop this confusion because the full Manifold can in fact be covered with a single copy of R^(N) so much of what we learn in differential geometry is redundant when applied to Euclidean space (or superfluous).  
When dealing with an arbitrary manifold, M, one maps open sets in R^(N) into open sets in M and then uses the standard coordinate in R^(N) to label points in M.  This labeling does not have the structure of a vector space.  Also, one cannot in general cover M with a single open set of R^(N).  While {r_i} is not a vector in M, differential displacements {dr_i} form a vector space.  They obey the axioms of a vector space.  Now, to make a differential displacement one needs a base point.  On a Manifold, M, we build this construction at each point of M.  The space made of all differential displacements at the point p is called the dual tangent space at p, T*(M, p) is one notation for this.  I say dual because the vector space at p is defined by the gradient of a scalar function on M evaluated at p, {df/dx_i}.  These quantities obey the general coordinate transform for the components of a vector in a vector space.  The tangent space of M at p is denoted T(M, p).   Using the standard coordinates from the mappings of open sets of R^(N) into M creates what we can a coordinate basis for T and T*.  We have a complete set of linearly independent elements generated by the {dx_i} and the {df/dx_i}.  Again, each point of M has its own tangent space.  
One can only do vector algebra operations on vectors in the spaces at the same point.  There is a one to one correspondence between elements of T and T* at each point of M.  The linear map between the two spaces is defined in terms of the metric tensor associated with M, g_ij.  With this one can define inner products (dot product generalization) as g_ij * V^i * U^j, sum over i and j, between two elements of T(M, p) or T*(M, p).  Typically, once we make the distinction between these two spaces we index elements of of one space with a lower index (subscript) and the other with an upper index (superscript).  This is equivalent to taking v_i * U^j, v_i = g_ik * V^k.  In other words we can project an element of T onto an element of T*.  If the underlying manifold is "flat", i.e. Euclidean with the typical metric = diag(1, 1, ..., 1) then T and T* are identical copies of each other, one cannot really decipher the difference between the two except as a formality.  
Now, in your case you state that the {q} are coordinates in the constraint manifold (or constrained manifold).  I can relate this to looking a systems in mechanics and mechanical engineering like a linkage, where there are objects connected by frictionless ideal joints that act as holonomic constraints.  As a simple example, a spherical pendulum would have a coordinate in R^(3), {x, y, z} but is constrained by r = constant.  Here the {q} would be {theta, phi}.  Your dq would be {d(theta), d(phi)}.  The Cartesian coordinates {x, y, z} are scalar functions of theta and phi.  Hence don't think of them as vectors in R^3 but as as set of scalar functions on S^2 (the 2-sphere).  Hence the derivatives, dx/d(theta) etc are elements of T(S^2, p), and dq are elements of T*(S^2, p).  Note that the indices line up properly, you are summing over the index of the {q} coordinates, for each r_i.  the d(r_i) on the l.h.s is also a dual vector.  Think of the dq as a basis for T*(S^2, p) and the partials of {r_i} as coefficients of d(r_i) in that basis.   
