Intersection between line and plane in R4-RN

[LMHmedchem]

Hello,

I believe that this is the correct forum for this post but please let me know if otherwise.

I have the following data in R4,

id     x              y             z            a
B     -0.0700177      0.382249     -0.289338     0.00349957
A      1.13388        1.77602      -0.679365    -0.729256
D      3.97525       -0.0977321    -0.319119     0.864731
E     -2.78909        0.10647      -0.206478    -0.0550962
C    -35.0653177    -39.836051     23.355262     1.57447957

This data is roughly in the form of, though this is in R3,

View attachment 8272

I need to find the coordinates for point F that is on line BC and also in the plane of A,D,E.

Generally, I guess that I would find the equation of the plane and the equation of the line and then the simultaneous solution for both equations. To find the equation of a plane, I would find the normal to the plane by taking the cross product of vectors AD and AE. Since the cross product is only defined in R3 and R7, I can't do this for R4, or any other dimensions.

I have vectors from point B to all other points and also vectors from point A to points D and E but I'm not sure how any of that helps. I assume that you can construct the equation for a hyperplane using a matrix operation, but I that's about as far as I have got with this.

I will code up the solution in c++ when I have a working method so I could use the Eigen library to to matrix operations if that is necessary.

Suggestions would be appreciated.

LMHmedchem

 

[I like Serena]

Hi LMHmedchem,

I'm afraid that a plane determined by 3 points is not a hyperplane in 4D. We need 4 points to define a hyperplane.

Consider that if we have two 2D planes in 4D, that this is similar to two lines in 2D.
Two lines in 2D can be either parallel, intersect in 1 point, or coincide.
Similarly two 2D planes in 4D can either be parallel, intersect in 1 point, intersect in a line, or coincide.

In this case you are intersecting a line with a 2D plane in 4D.
In other words, it is likely (I didn't check yet) that there is no intersection at all, just like a point in 2D is unlikely to be part of a line.

Can you clarify if we are supposed to have a hyperplane after all?
Or do you otherwise have a reason to expect that there will be a single point of intersection?

Oh, and I moved your thread to Linear Algebra.

 

[LMHmedchem]

Oh, and I moved your thread to Linear Algebra.

Thanks, I'm often a bit unsure about the best place to post.

I'm afraid that a plane determined by 3 points is not a hyperplane in 4D. We need 4 points to define a hyperplane.

Can you clarify if we are supposed to have a hyperplane after all?
Or do you otherwise have a reason to expect that there will be a single point of intersection?

I am not entirely sure about this. In my understanding, providing the space is orthogonal (which it is in this case) any three points can be considered to lie in a Euclidean plane, regardless of the dimensionality of the point coordinates. I would define this plane as two vectors with their tail at one of the three points such as vAD and vAE in my example above. I called it a hyperplane because the dimensions of the points are R4.

View attachment 8273

You can calculate the distance between A and D as the magnitude of the vector vAD. You can take the dot product between vectors vAD and vAE and thus find the angle. All of these calculations work for R4 or higher as well as for R2 and R3. You cannot, of course, find the cross product after R3, but you can do other operations in the plane such as finding vector projections and rejections. I have seen here and there other methods for calculating the equation of a plane that do not involve plane normals.

Based on this, I was thinking that you should be able to find the intersection between the plane vAD,vAE and a line vBC as in my example above. I am assuming that vBC is not parallel to the plane or in the plane. I would think that there would have to be a single point of intersection if neither of those are true.

Am I thinking about this incorrectly?

LMHmedchem

 

[I like Serena]

I called it a hyperplane because the dimensions of the points are R4.

That is not a hyperplane.
A hyperplane is a space that has one dimension less than the space it is contained in.
So in 4 dimensions a hyperplane is a 3-dimensional space.
It means that ADE is not a hyperplane but 'just' a plane.Anyway, let's see where we can get.

Let:
$$\mathbf v_1 = \mathbf D - \mathbf A\\ \mathbf v_2 = \mathbf E - \mathbf A\\ \mathbf v = \mathbf C - \mathbf B$$
and let $\boldsymbol \pi$ be the intersection point. That is, the projection of $\mathbf B$ in the direction of $\mathbf v$ onto the plane.

Then $\boldsymbol\pi$ must be $\mathbf A$ plus a linear combination of $\mathbf v_1$ and $\mathbf v_2$.
That is:
$$\boldsymbol \pi = \mathbf A + \lambda \mathbf v_1 + \mu \mathbf v_2 \tag 1$$
And we must also have that $\boldsymbol \pi$ is on the line BC.
That is:
$$\boldsymbol \pi = \mathbf B - \xi \mathbf v \tag 2$$
In other words:
$$\boldsymbol \pi = \mathbf A + \lambda \mathbf v_1 + \mu \mathbf v_2 =\mathbf B - \xi \mathbf v \quad\Rightarrow\quad
\lambda \mathbf v_1 + \mu \mathbf v_2 + \xi \mathbf v = \mathbf B - \mathbf A \quad\Rightarrow\quad
\begin{bmatrix}\mathbf v_1 & \mathbf v_2 & \mathbf v\end{bmatrix}\begin{bmatrix}\lambda \\ \mu \\ \xi\end{bmatrix} = \mathbf B - \mathbf A \tag 3
$$

These are 4 equations with 3 unknowns, meaning the system is overdetermined.
This is expected and means that (3) will generally not have a solution.
The best we can do is to find the point on the line BC that is closest to the plane ADE.

Let $M = \begin{bmatrix}\mathbf v_1 & \mathbf v_2 & \mathbf v\end{bmatrix}$.
Then the easiest way to do that, is to apply the so called Moore-Penrose pseudo inverse of $M$ to $\mathbf B - \mathbf A$.
That is:
$$\begin{bmatrix}\lambda \\ \mu \\ \xi\end{bmatrix} \approx M^+ (\mathbf B - \mathbf A) \tag 4$$
where $M^+$ is the pseudo inverse, and the approximation is such that we find the closest point to the plane.

Btw, instead of the pseudo inverse we can also complete the matrix to a 4x4 matrix by adding the vector that is perpendicular to the other 3 vectors, and take its inverse.
That's not particularly hard. I'm just keeping the formulas to a minimum for now.

We can calculate this for instance with Octave online:

B=[ -0.0700177   0.382249  -0.289338  0.00349957 ];
A=[  1.13388     1.77602   -0.679365 -0.729256   ];
D=[  3.97525    -0.0977321 -0.319119  0.864731   ];
E=[ -2.78909     0.10647   -0.206478 -0.0550962  ];
C=[-35.0653177 -39.836051  23.355262  1.57447957 ];
v1=D-A
v2=E-A
v=C-B
M=[v1' v2' v']
parm = M \ (B-A)'
lambda=parm(1)
mu=parm(2)
xi=parm(3)
pi=A + lambda * v1 + mu * v2
dist=norm(pi - (B - xi * v))

The result is:

v1 =
   2.84137  -1.87375   0.36025   1.59399
v2 =
  -3.92297  -1.66955   0.47289   0.67416
v =
  -34.9953  -40.2183   23.6446    1.5710
M =
    2.84137   -3.92297  -34.99530
   -1.87375   -1.66955  -40.21830
    0.36025    0.47289   23.64460
    1.59399    0.67416    1.57098
parm =
   0.2594223
   0.4660341
   0.0032226
lambda =  0.25942
mu =  0.46603
xi =  0.0032226
pi =
   0.0427568   0.5118597  -0.3655277  -0.0015588
dist =    9.0684e-06

In this particular case the distance of the closest point to the plane is pretty small (9.0684e-06).
But then, $\mathbf B$ was already pretty close to the plane ADE, which we can see since $\xi$ is pretty small too (0.0032226).

 

[LMHmedchem]

That is not a hyperplane.
A hyperplane is a space that has one dimension less than the space it is contained in.
So in 4 dimensions a hyperplane is a 3-dimensional space.
It means that ADE is not a hyperplane but 'just' a plane.

Thanks, my understanding always gets a bit fuzzy after R3.

I have a few clarification questions if you don't mind.

Question 1:

In the expression for the point we are looking for,
$$ \pi = \mathbf B - \xi \mathbf v $$
$$ \xi \text{ = some linear multiplier of } \vec {v} \text { that extends the line CB in the direction -B and ends in the plane ADE}$$
$$ \text {since } \xi \text { prescribes the distance BC must be extended, is this the magnitude of the vector } \vec {B\pi}$$
$$ \text {I guess I am asking if you knew the coordinates of } \pi \text{ could you back-calculate } \xi $$
$$ \text { in the manner you would calculate any other vector magnitude?}$$

Question 2:

These are 4 equations with 3 unknowns, meaning the system is overdetermined.
This is expected and means that (3) will generally not have a solution.
The best we can do is to find the point on the line BC is that closest to the plane ADE.

Does "generally not have a solution" mean,

a) In most cases no point will exist in the prescribed space that is both in the plane ADE and on the line BC.
b) A point will exist that is both in the plane ADE and on the line BC but there is no way to find the coordinates for that point (all we can find is an estimate).

Question 3:

In this particular case the distance of the closest point to the plane is pretty small (9.0684e-06).

Is this the distance of the point pi to the plane ADE or the distance of point B to the point pi?

Question 4:

Does the above work for points with any number of coordinates? It appears from your Octave code that more coordinates would just mean a larger matrix but the same solution would apply. Is that correct?

Thanks for the help so far, I will code this up and see if I can get it working.

I have some thoughts on a different solution that might be able to find a point that is demonstrably on the line and in the plane, but that of course would presuppose that such a point exists in every case so I will wait for your response to question 2.

LMHmedchem

 

[I like Serena]

Let me answer question 3 first and disambiguate the point in the plane ADE from the point on the line BC.

Question 3:

Is this the distance of the point pi to the plane ADE or the distance of point B to the point pi?

There is ambiguity between the point in the plane ADE that is closest to the line, and the point on the line BC that is closest to the plane. That's because we started out assuming that they were the same point. They are not.

Let's disambiguate them. Let $\boldsymbol\pi$ be the closest point in the plane ADE, and let P be the closest point on the line BC.
So:
$$\boldsymbol\pi = \mathbf A + \lambda \mathbf v_1 + \mu \mathbf v_2 \tag{1a}$$
and:
$$\mathbf P = \mathbf B - \xi \mathbf v \tag{2a}$$
The distance that I was talking about before is the distance between $\boldsymbol \pi$ and $\mathbf P$.

Question 1:

In the expression for the point we are looking for,
$$ \pi = \mathbf B - \xi \mathbf v $$
$$ \xi \text{ = some linear multiplier of } \vec {v} \text { that extends the line CB in the direction -B and ends in the plane ADE}$$
$$ \text {since } \xi \text { prescribes the distance BC must be extended, is this the magnitude of the vector } \vec {B\pi}$$
$$ \text {I guess I am asking if you knew the coordinates of } \pi \text{ could you back-calculate } \xi $$
$$ \text { in the manner you would calculate any other vector magnitude?}$$

Now that we have disambiguated $\boldsymbol\pi$ and $\mathbf P$, we can back-calculate $\xi$ from $\mathbf P$.
However, as $\xi$ is defined now, it's not the length of $\overrightarrow{BP}$.
That is because $\mathbf v$ does not have length $1$. We have to compensate for the length of $\mathbf v$.
Instead we have:
$$|\xi| = \frac{\|\overrightarrow{BP}\|}{\|\mathbf v\|}$$
where $\|\cdot \|$ is the norm function that gives the length of a vector.

Question 2:Does "generally not have a solution" mean,

a) In most cases no point will exist in the prescribed space that is both in the plane ADE and on the line BC.
b) A point will exist that is both in the plane ADE and on the line BC but there is no way to find the coordinates for that point (all we can find is an estimate).

The answer is (a). Generally no such point exists.

To make this clear, suppose $\mathbf n$ is a vector that is perpendicular to $\mathbf v_1, \mathbf v_2$, and $\mathbf v$.
Then we can extend the plane ADE into an actual hyperplane by including $\mathbf n$ in addition to $\mathbf v_1, \mathbf v_2$.
The corresponding equation becomes:
$$\mathbf P = \mathbf A + \lambda \mathbf v_1 + \mu \mathbf v_2 + \sigma \mathbf n = \mathbf B - \xi \mathbf v \tag{3a}$$
Now we have 4 equations with 4 unknowns, which does have a unique solution. So the line BC has a single point of intersection with the new hyperplane.
The distance to the plane ADE is then the distance along the new vector $\mathbf n$. If we pick $\mathbf n$ such that its length is $1$, this distance is $|\sigma|$.
Generally this distance will not be zero.

To understand this from an $\mathbb R^3$ point of view, consider 2 lines in $\mathbb R^3$. They will generally not intersect in a point will they?
However, if we extend one of the lines into a plane along a perpendicular vector, we will find a single point of intersection of the remaining line with the new plane.
Unless of course the lines were parallel or coincided to begin with.

Question 4:

Does the above work for points with any number of coordinates? It appears from your Octave code that more coordinates would just mean a larger matrix but the same solution would apply. Is that correct?

That is correct. The method works with an arbitrary number of coordinates.

 

[LMHmedchem]

Thanks for the explanations. I am going to code this up and see what results I can get but I have a few more clarification questions.

I am not sure what the single quote means in the Octave code you posted, such as,

M=[v1' v2' v']
parm = M \ (B-A)'

This looks like M is defined as a matrix of the three vectors v1, v2, and v but I don't understand what the ' signifies in v', etc. M looks like the transpose of the of the three vectors. Is that what the ' signifies?

It looks like the Octave code you posted is for the pseudo inverse solution. Is that correct? Can you also post the solution for the second solution you mentioned where you add the vector n to equation 3 of your second post so that the system is fully determined? What would the method be for finding vector n that is orthagonal to v, v1 and v2?

LMHmedchem

 

[I like Serena]

I am not sure what the single quote means in the Octave code you posted, such as,

M=[v1' v2' v']
parm = M \ (B-A)'

This looks like M is defined as a matrix of the three vectors v1, v2, and v but I don't understand what the ' signifies in v', etc. M looks like the transpose of the of the three vectors. Is that what the ' signifies?

Indeed. The ' is the transpose.
The vectors should really be column vectors. I made them row vectors only because it makes Octave's output more readible, and also because you more or less specified them as row vectors in post #1.

It looks like the Octave code you posted is for the pseudo inverse solution. Is that correct?

Yes. The solution of $M\mathbf x = \mathbf b$ in Octave is [M]x = M \ b[/M]. That is, we 'left divide' $\mathbf b$ by $M$.
If $M$ is a square invertible matrix, it is equivalent to $\mathbf x = M^{-1}\mathbf b$.
And if $M$ is either not square or not invertible, it is equivalent to $\mathbf x = M^{+}\mathbf b$.

Can you also post the solution for the second solution you mentioned where you add the vector n to equation 3 of your second post so that the system is fully determined? What would the method be for finding vector n that is orthagonal to v, v1 and v2?

To find a vector $\mathbf n$ that is perpendicular to $\mathbf v_1,\mathbf v_2$, and $\mathbf v$, we need to solve $\mathbf v_1 \cdot \mathbf n = 0, \mathbf v_2 \cdot \mathbf n = 0, \mathbf v \cdot \mathbf n = 0$.
In matrix notation:
$$\begin{cases}\mathbf v_1 \cdot \mathbf n = 0 \\ \mathbf v_2 \cdot \mathbf n = 0 \\ \mathbf v \cdot \mathbf n = 0 \end{cases} \quad\Rightarrow\quad
\begin{bmatrix}\mathbf v_1 & \mathbf v_2 & \mathbf v\end{bmatrix}^T \mathbf n = \mathbf 0 \quad\Rightarrow\quad
M^T\mathbf n = \mathbf 0 \tag{4a}$$
It means that $\mathbf n$ represents the null space of $M^T$. We can solve it for instance with Gaussian elimination. In Octave it's simply [M]n = null(M')[/M]. It will also pick $\mathbf n$ such that its length is $1$, so that the $|\sigma|$ in (3a) will be the distance between $\boldsymbol\pi$ and $\mathbf P$.

We can now solve (3a), which was:
$$\mathbf P = \mathbf A + \lambda \mathbf v_1 + \mu \mathbf v_2 + \sigma \mathbf n = \mathbf B - \xi \mathbf v \tag{3a}$$

Let $N = \begin{bmatrix}\mathbf v_1 & \mathbf v_2 & \mathbf v & \mathbf n\end{bmatrix}$.
It follows that:
$$N\begin{bmatrix}\lambda \\ \mu \\ \xi \\ \sigma\end{bmatrix} = \mathbf B - \mathbf A \tag{5a}$$
which we can again solve with for instance Gaussian elimination.
Finally we can find $\boldsymbol\pi$, $\mathbf P$, and their distance $|\sigma|$, by substituting what we found in (1a) and (2a).

In Octave:

n = null(M')
N = [M n]
parm = N \ (B-A)';
lambda=parm(1)
mu=parm(2)
xi=parm(3)
sigma=parm(4)
P = B - xi * v
pi = A + lambda * v1 + mu * v2
dist = abs(sigma)

with result:

n =
  -0.036705
   0.484310
   0.738375
   0.467867
N =
    2.841370   -3.922970  -34.995300   -0.036705
   -1.873752   -1.669550  -40.218300    0.484310
    0.360246    0.472887   23.644600    0.738375
    1.593987    0.674160    1.570980    0.467867
lambda =  0.25942
mu =  0.46603
xi =  0.0032226
sigma =   -9.0684e-06
P =
   0.0427571   0.5118553  -0.3655344  -0.0015630
pi =
   0.0427568   0.5118597  -0.3655277  -0.0015588
dist =    9.0684e-06

 

[LMHmedchem]

Hello Again,

I am in the process of writing this algorithm into c++ code to test and I have another question.

I have taken the transpose of the three vector matrix v1, v2, v (M=[v1' v2' v'])

B=[ -0.0700177   0.382249  -0.289338  0.00349957 ];
A=[  1.13388     1.77602   -0.679365 -0.729256   ];
D=[  3.97525    -0.0977321 -0.319119  0.864731   ];
E=[ -2.78909     0.10647   -0.206478 -0.0550962  ];
C=[-35.0653177 -39.836051  23.355262  1.57447957 ];

v1 = D-A
v2 = E-A
v  = C-B

v1 =   2.84137,  -1.87375,  0.36025,  1.59399
v2 =  -3.92297,  -1.66955,  0.47289,  0.67416
v  = -34.9953,  -40.2183,  23.6446,   1.5710

M = 2.84137,  -3.92297,  -34.99530
   -1.87375,  -1.66955,  -40.21830
    0.36025,   0.47289,   23.64460
    1.59399,   0.67416,    1.57098

The next step in the algorithm is to evaluate parm = M \ (B-A)'. If v3 = (B-A),

v3 = -1.20390, -1.39377, 0.39003, 0.73276

This expression divides the transpose M by v3' (transpose of B-A). Vector v3 has 4 elements, so its transpose would be a matrix with 1 column and 4 rows. This makes the expression parm = M \ (B-A)'

\[
param =
\left[ {\begin{array}{ccccc}
2.84137 & -3.92297 & -34.99530\\
-1.87375 & -1.66955 & -40.21830\\
0.36025 & 0.47289 & 23.64460\\
1.59399 & 0.67416 & 1.57098\\
\end{array} } \right]
{\LARGE /}
\left[ {\begin{array}{c}
-1.20390\\
-1.39377\\
0.39003\\
0.73276\\
\end{array} } \right]
\]

or,

\[
param =
\left[ {\begin{array}{ccccc}
2.84137 & -3.92297 & -34.99530\\
-1.87375 & -1.66955 & -40.21830\\
0.36025 & 0.47289 & 23.64460\\
1.59399 & 0.67416 & 1.57098\\
\end{array} } \right]
^{\LARGE *}
\left[ {\begin{array}{c}
-1.20390\\
-1.39377\\
0.39003\\
0.73276\\
\end{array} } \right]
^{\LARGE -1}
\]

Since dividing a numerator matrix by a denominator matrix is the equivalent of multiplying the numerator matrix by the inverse of the denominator matrix, it is necessary to take the inverse of the denominator matrix as a prerequisite step. I thought that it is only possible to take the inverse of a square matrix, so this is not adding up.

the following expression,

M = [ 2.84137,  -3.92297,  -34.99530;
     -1.87375,  -1.66955,  -40.21830;
      0.36025,   0.47289,   23.64460;
      1.59399,   0.67416,    1.57098 ]

M2 = [ -1.20390;
       -1.39377;
        0.39003;
        0.73276 ]

param = M/M2

results in the expected error at octave online, nonconformant arguments (op1 is 4x3, op2 is 4x1), so I obviously don't have this right.

Am I supposed to be dividing M by a vector instead of a 1x4? If so, what is the meaning of (B-A)'?

LMHmedchem

 

[I like Serena]

This makes the expression parm = M \ (B-A)'
\[
param =
\left[ {\begin{array}{ccccc}
2.84137 & -3.92297 & -34.99530\\
-1.87375 & -1.66955 & -40.21830\\
0.36025 & 0.47289 & 23.64460\\
1.59399 & 0.67416 & 1.57098\\
\end{array} } \right]
{\LARGE /}
\left[ {\begin{array}{c}
-1.20390\\
-1.39377\\
0.39003\\
0.73276\\
\end{array} } \right]
\]
or,
\[
param =
\left[ {\begin{array}{ccccc}
2.84137 & -3.92297 & -34.99530\\
-1.87375 & -1.66955 & -40.21830\\
0.36025 & 0.47289 & 23.64460\\
1.59399 & 0.67416 & 1.57098\\
\end{array} } \right]
^{\LARGE *}
\left[ {\begin{array}{c}
-1.20390\\
-1.39377\\
0.39003\\
0.73276\\
\end{array} } \right]
^{\LARGE -1}
\]

It's the other way around.
We're supposed to divide v3' on the left by M.
That would be $M^{-1}v_3^T$, except that obviously the inverse of $M$ does not exist because $M$ is not square.
The distinction comes from the fact that matrix multiplications are not commutative. That is, we cannot swap the arguments without affecting the result.
In my previous post I explained how we can construct the matrix $N$ such that we can evaluate $N \backslash v_3' = N^{-1} v_3^T$. And since $N$ is a square matrix, we can invert it.

the following expression,

M = [ 2.84137,  -3.92297,  -34.99530;
     -1.87375,  -1.66955,  -40.21830;
      0.36025,   0.47289,   23.64460;
      1.59399,   0.67416,    1.57098 ]

M2 = [ -1.20390;
       -1.39377;
        0.39003;
        0.73276 ]

param = M/M2

results in the expected error at octave online, nonconformant arguments (op1 is 4x3, op2 is 4x1), so I obviously don't have this right.

Am I supposed to be dividing M by a vector instead of a 1x4? If so, what is the meaning of (B-A)'?

It should be [M]param = M \ M2[/M], meaning we divide on the left instead of on the right.
Then the arguments will be conformant.

 

[I like Serena]

For the record, instead of [M]M \ v3'[/M] we can also write [M]pinv(M) * v3'[/M] in Octave.
It will yield the same result, since [M]pinv(M)[/M] calculates the pseudo inverse (or the real inverse if the matrix is actually invertible).
However, [M]M \ v3'[/M] is preferred since it allows Octave to use an optimized algorithm that does not necessarily calculate the pseudo inverse. Then again, for a 'small' matrix like the one we have here it won't matter.

 

[LMHmedchem]

For the record, instead of [M]M \ v3'[/M] we can also write [M]pinv(M) * v3'[/M] in Octave.
It will yield the same result, since [M]pinv(M)[/M] calculates the pseudo inverse (or the real inverse if the matrix is actually invertible).
However, [M]M \ v3'[/M] is preferred since it allows Octave to use an optimized algorithm that does not necessarily calculate the pseudo inverse. Then again, for a 'small' matrix like the one we have here it won't matter.

So Octave will automatically evaluate [M]M \ v3'[/M] with an inverse, pseudo inverse, or possibly some other appropriate solver, depending on what is efficient/possible based on the properties of M? Does Octave have a verbose output format that might explain more completely what it has done for a given problem?

I will work on the second solution you posted which I believe is using Gaussian elimination to find the vector n that is perpendicular to v1,v2, and v and then solving the system,

$$P=A+λv1+μv2+σn=B−ξv$$

I forgot to ask if this approach, like the moore-penrose solution, is the same no matter how many coordinate elements are in vectors v, v1, and v2. If each of these vectors has 20 elements instead of 4, do you still just need to add vector n to solve the system? I have interesting data with 50-200 elements per vector, so it is important that any approach I implement can at least possibly work in large dimensions.

LMHmedchem

 

[I like Serena]

So Octave will automatically evaluate [M]M \ v3'[/M] with an inverse, pseudo inverse, or possibly some other appropriate solver, depending on what is efficient/possible based on the properties of M? Does Octave have a verbose output format that might explain more completely what it has done for a given problem?

I will work on the second solution you posted which I believe is using Gaussian elimination to find the vector n that is perpendicular to v1,v2, and v and then solving the system,

$$P=A+λv1+μv2+σn=B−ξv$$

I looked it up, and apparently Octave will solve [M]M \ v3'[/M] with a least squares solution.
That is equivalent to solving:
$$(M^TM)\begin{bmatrix}\lambda \\ \mu \\ \xi \end{bmatrix} = M^T \mathbf v_3 \quad\Rightarrow\quad
\begin{bmatrix}\lambda \\ \mu \\ \xi \end{bmatrix} = (M^TM)^{-1}M^T \mathbf v_3
$$
Since $M^TM$ is a square matrix, it is invertible.

In other words, there is no need to find those perpendicular vectors, nor to find $N$ after all.

I forgot to ask if this approach, like the moore-penrose solution, is the same no matter how many coordinate elements are in vectors v, v1, and v2. If each of these vectors has 20 elements instead of 4, do you still just need to add vector n to solve the system? I have interesting data with 50-200 elements per vector, so it is important that any approach I implement can at least possibly work in large dimensions.

There would be more than 1 vector n.
So we would need to find $\mathbf n_1, ..., \mathbf n_{17}$.
Anyway, as I just wrote, that is not necessary after all.

 

[LMHmedchem]

That is equivalent to solving:
$$(M^TM)\begin{bmatrix}\lambda \\ \mu \\ \xi \end{bmatrix} = M^T \mathbf v_3 \quad\Rightarrow\quad
\begin{bmatrix}\lambda \\ \mu \\ \xi \end{bmatrix} = (M^TM)^{-1}M^T \mathbf v_3
$$
Since $M^TM$ is a square matrix, it is invertible.

I have been trying to code this up with the Eigen library but I seem to be misunderstanding the order of operations.

This expression,
$$ (M^TM)^{-1}M^T \mathbf v_3% $$
seems to involve either the product,
$$ ( (M^TM)^{-1}*M^T )* \mathbf v_3% $$
or the product,
$$ (M^TM)^{-1}*(M^T* \mathbf v_3)% $$

Since M^T is a 4x3 matrix, and v₃ is a 4x1 matrix, you cannot calculate the product,
$$ M^T* \mathbf v_3% $$

The (M^T*M)^-1*M^T matrix product is also a 4x3 matrix, so you cannot calculate that result product matrix by v₃ either.

I am guessing that I just have something out of order here.

My current data is,

matrix MT (transpose of M)
 2.8413699999999999513   -3.9229699999999998461   -34.995300000000000296
-1.8737500000000000266   -1.6695500000000000895   -40.218299999999999272
 0.36024600000000001065   0.47288700000000000179   23.644600000000000506
 1.593990000000000018     0.67415999999999998149    1.5709800000000000431

matrix (MT*M) ^-1 (inverse of M square matrix)
 -77717509171.034744263    1025437670022.9156494    1563374634806.2131348     990620125940.82434082
1025437670022.7844238    -13530058108103.224609   -20627825826984.5625     -13070660712841.994141
1563374634805.998291     -20627825826984.355469   -31449029630816.34375    -19927431979515.078125
 990620125940.70947266   -13070660712842.121094   -19927431979515.46875    -12626861622117.417969

matrix version of v3
 1.2038976999999999595
 1.3937710000000000932
-0.39002700000000001257
-0.73275557000000002272

What do I have wrong here?

LMHmedchem

 

[I like Serena]

I have been trying to code this up with the Eigen library but I seem to be misunderstanding the order of operations.

This expression,
$$ (M^TM)^{-1}M^T \mathbf v_3% $$
seems to involve either the product,
$$ ( (M^TM)^{-1}*M^T )* \mathbf v_3% $$
or the product,
$$ (M^TM)^{-1}*(M^T* \mathbf v_3)% $$

It does not matter which product is evaluated first. You'll see.

Since M^T is a 4x3 matrix, and v₃ is a 4x1 matrix, you cannot calculate the product,
$$ M^T* \mathbf v_3% $$
The (M^T*M)^-1*M^T matrix product is also a 4x3 matrix, so you cannot calculate that result product matrix by v₃ either.

$M$ is supposed to be a 4x3 matrix, so that $M^T$ is a 3x4 matrix.
Consequently $M^TM$ is a 3x3 matrix, and $(M^TM)^{-1}$ is a 3x3 matrix as well.

 

[LMHmedchem]

$M$ is supposed to be a 4x3 matrix, so that $M^T$ is a 3x4 matrix.
Consequently $M^TM$ is a 3x3 matrix, and $(M^TM)^{-1}$ is a 3x3 matrix as well.

I forgot to transpose the three vectors v1, v2, and v when I made M.

Thanks,

LMHmedchem

 

[LMHmedchem]

I do have this working now using the Eigen library and it is pretty fast. Let me know if anyone would like the c++ code posted.

I think I am still unsure about calculating the distance between P and pi, P and B, and pi and B. Normally I would just take the vector between each pair of points and then the magnitude of the vector would be the distance. Is this correct or do I need to do more than that?

LMHmedchem

 

[I like Serena]

I think I am still unsure about calculating the distance between P and pi, P and B, and pi and B. Normally I would just take the vector between each pair of points and then the magnitude of the vector would be the distance. Is this correct or do I need to do more than that?

That's indeed the way to find the distance between 2 points.