Re: Feature Request: Left (and Right) inverse of the Cross Product

On Mon, 2 Jun 2025 at 14:26, Leon Deligny via NumPy-Discussion
<[email protected]> wrote:
>
> ### Proposed new feature or change:
>
> Motivations: This is specific to 3D vector algebra. In Fluid Dynamics, we have access to the moment of a force at a specific point (M_P = OP \cross F). This calculation is crucial when determining the center of pressure (CoP), a pivotal concept for understanding the distribution of forces on an object submerged in a fluid. To accurately pinpoint the CoP, we often need to "reverse" this process, effectively requiring an inverse functionality for the cross product.
>
> How to compute the right-inverse: Given two real numbers a and c we want to find b such that c = a \cross b . If we write this equation in matrix format then we need to define:
>
> ```latex
> A = \begin{pmatrix}
> 0 & -a_3 & a_2 \\
> a_3 & 0 & -a_1 \\
> -a_2 & a_1 & 0
> \end{pmatrix}
> ```
>
> to get $c = A \cdot b$ (where $\cdot$ is the matrix multiplication). In real case scenarios there does not exist a vector b such that $c = A \cdot b$. So we can always find a vector b such that $|c - A \cdot b|_2^2$ is minimal (i.e. b is the best approximation $c \approx A \cdot b$).

The pseudo inverse of a skew symmetric 3x3 matrix like this is
proportional to itself. Here is a demonstration using sympy:

In [129]: x, y, z = symbols('x, y, z', real=True)

In [130]: M = Matrix([[0, -z, y], [z, 0, -x], [-y, x, 0]])

In [131]: print(M)
Matrix([[0, -z, y], [z, 0, -x], [-y, x, 0]])

In [132]: print(simplify(M.pinv()*(-x**2 - y**2 - z**2)))
Matrix([[0, -z, y], [z, 0, -x], [-y, x, 0]])

Given the matrix A then your least squares solution computed using numpy is:

In [133]: (A @ c) / ((A**2).sum() / -2)
Out[133]: array([-0.93244529,  0.74775365,  1.09117371])

That matches what Robert showed with lstsq:

In [134]: np.linalg.lstsq(A, c)[0]
Out[134]: array([-0.93244529,  0.74775365,  1.09117371])

Note that ((A**2).sum() / -2) is just -|a|^2 so if a is a unit vector
then this part can be simplified.

In vectors the minimum norm solution for b in a x b = c is just bhat =
(c x a) / |a|^2 so the cross product computes its own "inverse":

In [135]: np.cross(c, a) / (a**2).sum()
Out[135]: array([-0.93244529,  0.74775365,  1.09117371])

--
Oscar
_______________________________________________
NumPy-Discussion mailing list -- [email protected]
To unsubscribe send an email to [email protected]
https://mail.python.org/mailman3//lists/numpy-discussion.python.org
Member address: [email protected]