24 replies to this topic

[Mihail121]

Hi, consider the following problem: given a vector v=(x,y,z) located in a hyperplane. Given two reference vectors from that plane v1 and v2. Find a and b such that v = a*v1 + b*v2. How would you do that? a, b exist of course since any vector of a plane has a linear combination of 2 non-colinear vectors from the same plane. My maths are also not rusty and I generally know how to solve the system above, i.e.

x = a*x1 + b*x2
y = a*y1 + b*y2
z = a*z1 + b*z2

But what bothers me is this: to find a and b with the above system you need at most 2 of the 3 equations. But which two? Taking an equation out means you're projecting the 3D vectors on a 2D plane, which might give a headache, if after the projection the vectors become colinear, which is very much possible. So before doing the above calculation one has to check on which plane the reference vectors are not colinear and then choose the right 2 equations (cooridnates). There's just gotta be an easier and not that ugly way to do this I think, does anyone have an idea? Thank you!

[macnihilist]

I'm not entirely sure and I have not tried it, but the pseudoinverse of the 3x2 matrix

[v1 v2]

should project a vector v into the coordinate system [v1,v2] 'as good as possible' and give you a and b.

For real-valued matrices the Moore-Penrose-Pseudoinverse is (M^T * M)^-1 * M^T. (If the columns of M are linearly independent.)

[Mihail121]

Ah, interesting! Never knew about that property of the pseudo-inverse. I would definetely try it. By the way, before I sit down tonight and work out the maths, can you recommend me a good reference, describing these features of the pseudo-inverse?

edit: Are all types of pseudo-inverse matrices suitable or just the Moore-Penrose?

[macnihilist]

I just tried it in MATLAB with a few simple cases and it seems to work.

I don't have an online reference ready, but I think about it that way:
The pseudoinverse is a way to 'solve' the overdetermined least squares problem

min norm(A*x - d)

In your case

x = [a;b]
A = [v1,v2]
d = v

One solution is then

x = pinv(A)*d

But again: Be careful, this is based on the fading memories of a math lecture I had years ago...

[EDIT: of course it is pinv(A)*d, not pinv(A)*b]

[Mihail121]

Ich danke dir herzlichst! :)

[Nick]

There's a far more intuitive way to solve this.

Construct a vector v3, perpendicular to v1 and v2. This can be done using a cross product. Now you know that v has to be a linear combination of v1, v2, and v3, with the coefficient c for v3 in theory being 0:

x = a*x1 + b*x2 + c*x3
y = a*y1 + b*y2 + c*y3
z = a*z1 + b*z2 + c*z3

This is an ordinary system of three linear equations with three variables. Even if c isn't exactly 0 due to numerical imprecision, you know a and b are the coefficients of the projection of v onto the plane.

[.oisyn]

You can also approach this problem geometrically. In order to get to v, you start at the origin, and then go in the direction of v1 for a certain amount (a), and from thereon in direction of v2 which points directly at v in another amount (. The point where you change directions is the intersection p of the line with direction v1 through the origin, and that of the line with direction v2 through v. p will then be a*v1, and (v-p) will be b*v2. Calculating a and b should be trivial.

[Nick]

.oisyn said:

The point where you change directions is the intersection p of the line with direction v1 through the origin, and that of the line with direction v2 through v.

Problem is you can't easily intersect two lines in 3D. You could however intersect a line and a plane, with the plane making use of v3, the cross product of v1 and v2. That makes things no easier than my solution though. The intersection point has to be computed using a linear system of equations, after which you still need to determine a and b.

[Mihail121]

@.oisyn:

I need to do this very quickly and many many BILLION times per second... The geometrical interpretation does not help in that account, but it was the way I came into these calculations in the first place.

@.nick:

The extended system seems to give a solution in this particular case, yes, but seems a bit computationally intensive, I don't know, I have to benchmark it and I will. Also not sure if it's ill-conditioned and/or stable. The way I was about to do it, before you posted, was by extracting the least squares approximation of the overdetermined system (as macnihist suggested), but not using the pseudoinverse (best solution, but computationally intensive). Instead I'm about to calculate it using QR decomposition. The decomposition itself is performed using Givens rotations (see the great book "Matrix Computations" by Gehe H. Golub et. al.) Let's see if it works.

[coelurus]

If you can cache the lengths of v1 and v2, or even better always have them normalised, cheap vector projection would do the trick. If I missed any mention of that property, well, don't bother telling me :p

[.oisyn]

Nick said:

Problem is you can't easily intersect two lines in 3D. You could however intersect a line and a plane, with the plane making use of v3, the cross product of v1 and v2.

Yes, a line line intersection would obviously imply that one of the lines is converted into a plane.

Quote

That makes things no easier than my solution though.

I never said it was easier, it's just another way of looking at the problem . It might be easier for the intuition though, as with most problems solved geometrically rather than algebraically.

Mihail121 said:

@.oisyn:

I need to do this very quickly and many many BILLION times per second... The geometrical interpretation does not help in that account, but it was the way I came into these calculations in the first place.

Ah I see. How volatile are v, v1 and v2?

[Mihail121]

coelurus said:

If you can cache the lengths of v1 and v2, or even better always have them normalised, cheap vector projection would do the trick. If I missed any mention of that property, well, don't bother telling me :p

HANS, YOU'RE BACK!!!!!!!!!!!! I'm so glad to know you're alive :) Yes, the reference vectors (base of the plane) are normalized.

.oisyn said:

Ah I see. How volatile are v, v1 and v2?

I don't know yet, but my (crude) guess is they are VERY volatile. We shall see after some experiments.

[Nick]

Mihail121 said:

The extended system seems to give a solution in this particular case, yes, but seems a bit computationally intensive...

It's the least computationally intensive solution.

Quote

Also not sure if it's ill-conditioned and/or stable.

It's robust as long as v1 and v2 are not zero vectors and don't coincide (i.e. they form a well-defined plane).

Quote

Yes, the reference vectors (base of the plane) are normalized.

In that case the determinant is 1, which makes computing the 3x3 matrix inverse cheaper. Note that in fact you only need 2 / 3 of the matrix inverse, since you're not interested in the value of c.

[coelurus]

Hi Mikhail, I have been lurking on these forums every now and then, never dared to post anything :p

To project vector A onto vector B, you do "A' = (A.B) / (B.B) * B". If B is normalised, you have just "A' = A.B * B".

[Nick]

coelurus said:

To project vector A onto vector B, you do "A' = (A.B) / (B.B) * B". If B is normalised, you have just "A' = A.B * B".

That's a perpendicular projection. I don't think that's of use here.

[coelurus]

That was silly of me.
I have completely forgotten the instability behaviour with least squares fitting, would the normal solution really be a problem in this case?
QR decomposition definitely sounds overkill.

[macnihilist]

@Nick: Isn't the determinant only 1 if v1 and v2 are also perpendicular?

Anyway:
Solving Nick's system directly for a,b is probably the fastest (and best) solution, but after massaging the pseudoinverse stuff a bit further I think I have come up with a way that does not employ the 'full construction' of the pinv and offers a nice geometrical interpretation.
Of course this can also be seen as a solution to Nick's linear system, depends on how you approach it.

The rows of the pseudoinverse are two normals onto which you project v:

a = dot(n1,v)
b = dot(n2,v)

you have three constraints on each normal:

dot(n1,v1) = 1
dot(n1,v2) = 0
dot(n1,cross(v1,v2)) = 0 (distances out of the plane should not count)

instead of solving for n1 and n2, in this special case you can construct them directly:

c1 = cross(v1,v2);
n1 = cross(v2,c1);
n1 *= 1/dot(v1,n1);

c2 = cross(v2,v1); // or -c1
n2 = cross(v1,c2);
n2 *= 1/dot(v2,n2);

Then you have, like written above:

a = dot(n1,v)
b = dot(n2,v)


I tried it in MATLAB with 100000 random vectors and it worked. (Let's hope my program isn't flawed...)

[Mihail121]

coelurus said:

That was silly of me.
I have completely forgotten the instability behaviour with least squares fitting, would the normal solution really be a problem in this case?
QR decomposition definitely sounds overkill.

Yes, you were always silly :> Just don't get too smart around that particle accelerator :D Ok, back to the topic, the normal equation you mean? I never really considered that as an option, not a single man seems to be fond of that, because of extreme instability. QR decomp. is not really an overkill if fast Givens or, hmm, "modified" Givens is used.

Of course I don't have a problem to get the vectors normalized AND orthogonal and just transpose the coef. matrix to get the inverse, which will give me an answer instantly.

I was really hoping this thing has a small solution (1 mul, 2 adds or something), but that's life and it's complex.

@mac

Just tried it also, seems to work fine.

[Nick]

macnihilist said:

Isn't the determinant only 1 if v1 and v2 are also perpendicular?

Ah, right. The determinant of [v1 v2 (v1 x v2)] is (v1 x v2).(v1 x v2), which is the same as v3.v3. So it's not 1 but still a very simple computation.

[Nick]

Mihail121 said:

I was really hoping this thing has a small solution (1 mul, 2 adds or something), but that's life and it's complex.

I'm curious why you have to do this in the first place. If v is known to be in the plane, it was constructed as a linear combination of v1 and v2 in the first place. So instead of working in world xyz coordinates, why not work in the plane's coordinate system instead (ab coordinates)?

Can you reveal anything more about the actual application you need this for?