On Mon, 2 Jun 2025 at 14:26, Leon Deligny via NumPy-Discussion
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 -- numpy-discussion@python.org
To unsubscribe send an email to numpy-discussion-le...@python.org
https://mail.python.org/mailman3//lists/numpy-discussion.python.org
Member address: arch...@mail-archive.com
