Friday, March 19, 2021

The many faces of projection

What I want to talk about in fact is the Jacobian matrix of the constraint function. All these things that I want to mention are probably clear to people with strong backgrounds in linear algebra and differential geometry. But for me they were revelations of sort many of them.

All right, so maybe the first thing to start with is the fact the Jacobian matrix $\mathbf{J}$ of a function $g$ is the analog of the gradient operator or the normal to a surface. And this is why for many people $\mathbf{J}$ is synonymous to the constraint direction(s). Because the normal to the constraint manifold is transversal to the space of motion, hence along the constrained directions that can't move. Think of a pendulum: the constraint manifold is a circle and the normal to it is along the radius, which is exactly the pendulum's (constrained) direction.

But some math books also call $\mathbf{J}$ the differential (in coordinates). And this takes me to a confusion I always made in the beginning: for a single constraint, should $\mathbf{J}$ be a row or column vector? And it seemed other people weren't sure about it either in their notation. But I stuck with it being a row vector, as that was the majority view. And now it makes sense, because it is a covector (mainly because it has the components of a differential which is a 1-form, same thing). Hence, $\mathbf{J} \equiv dg$ and $\mathbf{J}^T \equiv \nabla g$.

I still don't fully grasp the concepts and that's why I don't want to talk formally about these things. But I do have to stress the fact that the configuration space of our mechanical system (where all the generalized coordinates reside) is not really $\mathbb{R}^n$, but one that has a metric (a curved space) given by the mass matrix $\mathbf{M}$ (i.e., the kinetic energy). It basically tells you that not all bodies are equal but weighed by their mass information. If this weren't the case, then vectors and covectors in $\mathbb{R}^n$ would be the same thing. But they're not! If for velocity (tangent) vectors $\mathbf{v}$ the norm is given by the well known kinetic energy ($\tfrac{1}{2}\mathbf{v}^T \mathbf{Mv}$), then for covectors the matrix is now the inverse of  $\mathbf{M}$. Think about momenta $\mathbf{p}$, the kinetic energy is the same, but now we have to divide by the mass rather than multiply ($\tfrac{1}{2}\mathbf{p}^T \mathbf{M}^{-1}\mathbf{p}$). And you have to multiply velocities by mass to get momenta $\mathbf{p} = \mathbf{Mv}$ and the other way around. Some other covectors are forces (as they work directly on the momenta) and surface normals; and they inhabit a place called the cotangent space.

Now let's fall back from geometry to linear algebra. Locally, $\mathbf{J}$ is a linearization of $g$ and maps the tangent space of the configuration space to the tangent space of the "constraint" space $\mathbb{R}^m$. So, it's a linear transform, thus a rectangular matrix. And it follows all the rules of linear algebra. Ideally, it's full rank and fat (rather than thin, i.e., $m <n$) in order for the system to be solvable. The image of  $\mathbf{J}$ is the whole of $\mathbb{R}^m$, which can be viewed as just a collection of numbers, just as for one constraint you get only one number. The meaning of these numbers comes from the analogy with the normal vector, in fact any vector for that matter. They are the magnitudes or components of the projection of a vector $\mathbf{v}$ on the "normal" $\mathbf{J}$, so $\mathbf{Jv}$ is just that. For one constraint it is a number, in general it is a $m$-dimensional vector because the "normal" $\mathbf{J}$ spans an $m$-dimensional space. You can think of taking each row of $\mathbf{J}$ and dotting it with $\mathbf{v}$ to get the projection length. And in fact, from linear algebra we know that the row space of $\mathbf{J}$ is formed by vectors $\mathbf{J}^T\mathbf{y}$, so we know that the rows of $\mathbf{J}$ span the "normal" or "constraint" space of dimension $m$. This space is orthogonal to the tangent space of the constraint manifold, which is precisely the null space of $\mathbf{J}$, i.e., vectors satisfying $\mathbf{Jv}=0$.

What does this tell us about $\mathbf{J}$? That it takes any vector $\mathbf{v}$ and converts it into its projection components. What about the reverse direction? It will only reconstruct the "normal" part of those vectors. So, you lose some information, but this is fine because you obtain the projected vector $\mathbf{J}^T\mathbf{Jv}$ (almost, stay tuned). If you think about it, this is what you do when you project on a unit vector $\mathbf{n}$: $\mathbf{nn}^T\mathbf{v}$. So you may say that $\mathbf{J}^T\mathbf{Jv}$ is a projection operator. Note that if $\mathbf{v}$ is already in the row space of $\mathbf{J}$ no information is lost and projection leaves it unchanged. Some people call $\mathbf{Jv}$ a "projection" and $\mathbf{J}^T\mathbf{v}$ an "un-projection", but this is not really true, as they do not satisfy the definition of a projection operator; but they do part of the job, so we may call these matrices something like "projectors" (when in fact they define the space we are projecting upon).

But this not the whole truth. In 3D you also have to divide by the norm squared $\mathbf{n}^T\mathbf{n}$ if we're dealing with a non-unit vector. And you have to do the same in $n$-dimensional space. Linear algebra and least square problems give you the formula $\mathbf{A}^{-1}\mathbf{Jv}$ for computing the actual projection lengths, where $\mathbf{A}=\mathbf{JJ}^T$. Here $\mathbf{A}$ plays the role of the norm squared of  $\mathbf{J}$ that you divide by. Finally, the projection that takes these lengths back to configuration space is $\mathbf{P}=\mathbf{JA}^{-1}\mathbf{J}^T$ and this is a genuine projection operator. Of course, if you're lucky and $\mathbf{A}$ is the identity then the previous paragraph still works, but that's rarely the case.

The projection $\mathbf{P}$ or rather its complement $\mathbf{I}-\mathbf{P}$ is usually called the "canonical" projection on the constraint manifold. And you may recognize the formula from the analytical solution of linearized velocity-level constrained dynamics. In practice you usually do not invert $\mathbf{A}$, but rather solve a linear system where the unknowns are the projection lengths. But $\mathbf{A}$ is not really the same, is it? 

So, what is the meaning behind $\mathbf{A}$ anyway? Some people call it the reduced or effective mass, others call it the Delassus operator, others compliance and so on. I think I hinted already at it: $\mathbf{A}$ determines the metric of the "constraint" space. In fact, it's the inverse $\mathbf{A}^{-1}$ that is the true effective mass, the one seen by the degrees of freedom along the constraint directions, and helps define the kinetic energy there. (I picture it as the transformation of the mass matrix $\mathbf{M}$ into "constraint" space). And $\mathbf{A}$ is actually the metric of the cotangent "constraint" space. This makes sense as $\mathbf{J}$ spans the cotangent space of the constraint manifold (in configuration space) and there the metric is $\mathbf{M}^{-1}$, so our "norm squared" becomes $\mathbf{A}=\mathbf{JM}^{-1}\mathbf{J}^T$. This is exactly the matrix that appears in constrained dynamics.

I want to end by noting the fun fact that by deriving a projection operator (pretty similar to what you would do in 3D), we actually solved for the Lagrange multipliers in constrained dynamics. In fact, projection is an optimization problem: find the shortest distance from a point to a surface. And the closest point happens to be on the perpendicular to the surface, one like our "normal" $\mathbf{J}$. Pretty much what we talked about. So, deriving the projector $\mathbf{P}$ amounts to solving the projection (or least squares) problem. Or it does at least for velocity vectors. The general case for projecting a point in configuration space on the constraint manifold is a bit trickier as it's a nonlinear problem, but it can be broken down into such linear steps. But that's another story...

To summarize, in order to solve the (linearized) constrained dynamics problem one simply needs to project velocities into the tangent space which is the null space of $\mathbf{J}$. This is done by projecting them onto the space spanned by the columns of $\mathbf{J}^T$ and then taking the orthogonal complement of that projection.

No comments:

Post a Comment