[Numpy-discussion] Feature Request: Left (and Right) inverse of the Cross Product
### 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$).
When minimizing $|c - A \cdot b|_2^2$ there exists an analytical solution found
with gaussian elimination procedure (we are working with 3x3 matrices and
3-dimensional vectors).
I did not find any litterature on this subject specifically. But this would be
a greatly appreciated feature.
___
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
[Numpy-discussion] Re: Feature Request: Left (and Right) inverse of the Cross Product
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
[Numpy-discussion] Re: Feature Request: Left (and Right) inverse of the Cross Product
On Mon, Jun 2, 2025 at 1:09 PM Oscar Benjamin via NumPy-Discussion < numpy-discussion@python.org> wrote: > 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]) > Delightful! -- Robert Kern ___ 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
[Numpy-discussion] Re: Feature Request: Left (and Right) inverse of the Cross Product
On Mon, Jun 2, 2025 at 9:30 AM Leon Deligny via NumPy-Discussion <
numpy-discussion@python.org> 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$).
>
> When minimizing $|c - A \cdot b|_2^2$ there exists an analytical solution
> found with gaussian elimination procedure (we are working with 3x3 matrices
> and 3-dimensional vectors).
>
> I did not find any litterature on this subject specifically. But this
> would be a greatly appreciated feature.
>
C.f. https://github.com/numpy/numpy/issues/29110
To be clear, I still don't think we're entertaining a feature request for
inclusion in numpy, but we on the mailing list can perhaps help you
implement this in your own code.
Note that because `A` is singular, there are an infinite number of `b`
values that will meet the criterion. The problem is not reducing the
residual (you can always get to 0, or near enough in floating point terms),
the problem is picking which of the infinite valid solutions you want.
`np.linalg.lstsq()` will give you the minimum-norm solution, which is
usually a pretty good option, though it may or may not suit your specific
application.
>>> import numpy as np
>>> rng = np.random.default_rng(0x7cadf78fc1ad960cdb91bb4745b84b28)
>>> av = rng.multivariate_normal(np.zeros(3), np.eye(3), size=10)
>>> bv = rng.multivariate_normal(np.zeros(3), np.eye(3), size=10)
>>> cv = np.cross(av, bv)
>>> cv
array([[ 0.78079691, -0.61265949, 0.43715586],
[ 0.0508124 , 2.3512873 , 0.17093971],
[ 0.39231732, 0.35281331, 0.90361588],
[-0.36610523, 0.24177299, 0.49421533],
[-0.51625536, 1.69312288, -0.9220552 ],
[-0.59298882, 0.60810492, 2.38083992],
[ 0.17308949, -0.59307124, -0.7205352 ],
[-0.17391941, -0.43438387, -0.451717 ],
[ 1.15385801, -0.80682863, -1.4479563 ],
[-1.77839153, -2.30249409, 0.05814433]])
# Fancy way to compute the A matrix corresponding to the operation
`np.cross(a, ...)
>>> A = np.cross(av[0], -np.eye(3))
>>> A @ bv[0]
array([ 0.78079691, -0.61265949, 0.43715586])
>>> cv[0]
array([ 0.78079691, -0.61265949, 0.43715586])
>>> for a, b in zip(av, bv):
... c = np.cross(a, b)
... A = np.cross(a, -np.eye(3))
... bhat = np.linalg.lstsq(A, c)[0]
... res = c - np.cross(a, bhat)
... print(res)
...
[-2.22044605e-16 0.e+00 -1.66533454e-16]
[7.56339436e-16 8.88178420e-16 2.77555756e-16]
[-5.55111512e-17 -2.77555756e-16 -2.22044605e-16]
[-1.66533454e-16 -5.55111512e-17 0.e+00]
[0.e+00 6.66133815e-16 3.33066907e-16]
[-1.11022302e-16 0.e+00 -4.44089210e-16]
[ 0.e+00 -1.11022302e-16 -2.22044605e-16]
[-1.11022302e-16 1.11022302e-16 -5.55111512e-17]
[ 0.e+00 0.e+00 -2.22044605e-16]
[2.22044605e-16 4.44089210e-16 0.e+00]
--
Robert Kern
___
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
