topic Re: How to calculate Vector Cross Product in Discussions

How to calculate Vector Cross Product

CaseyL — Mon, 25 Nov 2019 15:23:00 GMT

From the documentation it looks like the * operator will perform a cross product and ;* will perform a dot product. I want to use a cross product to determine the vector normal to a pair of known vectors like this:

p1 = [1, 2, 4];
p2 = [2, 1, 4];
p3 = [2, 2, 4];

v1 = p3 - p1;
v2 = p2 - p1;

xProd = v1*v2;

3 points define a plane, two vectors taken from those points do the same, the cross-product should give the normal vector of that plane. The problem seems to be that the dimensions of the two matrices being multiplied (v1 and v2) do not agree -- nRows(v1) should equal nCols(v2). I can implement this manually with a simplified formula for vector cross products but does anyone know how to format this differently to work with matrix multiplication?

Re: How to calculate Vector Cross Product

cwillden — Mon, 25 Nov 2019 16:58:35 GMT

If you're only ever dealing 3x1 (or 1x3) vectors, then I would think computational efficiency or elegance is not too important. Would this be sufficient?

cross_prod = function({v1, v2},
	matrix({v1[2]*v2[3] - v1[3]*v2[2], v1[1]*v2[3] - v1[3]*v2[1], v1[1]*v2[2] - v1[2]*v2[1]});
);

cross_prod(v1, v2);

Re: How to calculate Vector Cross Product

Craige_Hales — Mon, 25 Nov 2019 16:55:16 GMT

I hope someone comes up with a clever way to do this. The det() function is closely related, but I don't understand how to use det(3x3 matrix) to get back more than a scalar value. I think you could use det(2x2 sub-matrix) three times, but I think the straight-forward xyzzy approach is simpler and just as fast.

I use it for the surface normals in 3D scenes, there is some JSL in https://community.jmp.com/t5/Uncharted/Custom-Visualization/ba-p/191559 .

Re: How to calculate Vector Cross Product

CaseyL — Mon, 25 Nov 2019 18:29:47 GMT

Thanks cwillden, this is essentially what I have implemented. It does work for me since I am only concerned with 3D vectors, but I am supprised this isn't a built in function. I'm going to leave this open for a bit to see if anyone else knows of such a way.

Re: How to calculate Vector Cross Product

klk — Thu, 09 Mar 2023 19:24:34 GMT

Sorry for reviving this zombie thread, but there is a sign error in the second element here. I believe it should be

Matrix(
		{v1[2] * v2[3] - v1[3] * v2[2], v1[3] * v2[1] - v1[1] * v2[3], v1[1] * v2[2] - v1[2] * v2[1]}
)

Just want to leave this here for anyone else who comes along for a quick copy-paste.

Re: How to calculate Vector Cross Product

brady_brady — Tue, 14 Mar 2023 05:19:15 GMT

Well, the revived zombie caught me. Here is a way to do it via determinants as @Craige_Hales mentioned. It is short, but I don't know... it isn't very simple.

crossP = Function( {x, y, i = 1},
	Return( J( 3, 1, Det( (x || y)[Remove( [1, 2, 3], i++ ), 0] ) ) )
);

i.e., using the originally-supplied data:

Re: How to calculate Vector Cross Product

Craige_Hales — Tue, 14 Mar 2023 11:43:48 GMT

@brady_brady Cool! It needs a tweak if you are going to use it:

	i = 1;
	j = -1;
	J( 3, 1, (j *= -1) * Det( (v1 || v2)[Remove( [1, 2, 3], i++ ), 0] ) );

https://mathinsight.org/cross_product_formula shows how the middle element needs to be negative. Above, j will alternate sign (1, -1, 1), and i increases (1, 2, 3). The J(...) function has an implicit loop that evaluates the 3rd argument for each element.

The straight forward corrected formula from @klk is about 3 times faster and about N times easier to understand.

Test code.

v1 = [0, 0, 1];
v2 = [0, 1, 0];

start = HP Time();
For( t = 1, t < 1e6, t += 1,
	i = 1;
	j = -1;
	J( 3, 1, (j *= -1) * Det( (v1 || v2)[Remove( [1, 2, 3], i++ ), 0] ) );
);
stop = HP Time();
Show( (stop - start) / 1e6 );//9.3s


start = HP Time();
For( t = 1, t < 1e6, t += 1,
	Matrix( {v1[2] * v2[3] - v1[3] * v2[2], v1[3] * v2[1] - v1[1] * v2[3], v1[1] * v2[2] - v1[2] * v2[1]} )
);
stop = HP Time();
Show( (stop - start) / 1e6 );//3.2s


For( q = 1, q < 10000, q += 1,
	v1 = J( 3, 1, Random Uniform( -1, 1 ) );
	v2 = J( 3, 1, Random Uniform( -1, 1 ) );
	If(
		All(
			Round( Matrix( {v1[2] * v2[3] - v1[3] * v2[2], v1[3] * v2[1] - v1[1] * v2[3], v1[1] * v2[2] - v1[2] * v2[1]} ), 10 ) //
			== //
			(i = 1; j = -1 ; Round( J( 3, 1, (j *= -1) * Det( (v1 || v2)[Remove( [1, 2, 3], i++ ), 0] ) ), 10 ) ; ) //
		) == 0//
	, //
		Throw( Char( q ) )
	);
);

Curiously, rounding to 11 places occasionally fails to match results. Probably nearly collinear vectors.

Re: How to calculate Vector Cross Product

brady_brady — Tue, 14 Mar 2023 15:02:52 GMT

The inferior performance doesn't surprise... thanks for inserting the sign adjustment. I omitted the "sign checkerboard" inadvertently.

IMO this is one of those instances where "just because you can, doesn't mean you should" applies... {cool, fun, short} doesn't outweigh {understandable, maintainable, speedy}.