The intersection of M hyperplanes (Ndim)

Here is my understanding of what the theory behind what I am trying to do is, correct me where I am wrong, also please check my questions at the end:

The intersection of N hyperplanes (Ndim) is also known as the solution to "the system of N linear equations(=planes) with N variables". However, I am interested in the case when I have less equations. The result is not a single point in Ndim but a lower-dimensions plane. For example in 3D. The intersection of 3 planes (M=3, N=3), assuming any pair is not parallel, is a single 3D point and the intersection of 2 planes (M=2, N=3) is a line. In Ndim, the solution would be an Mdim hyperplane "living" in Ndim. Below, I call this a "hyperline".

An example of the above problem is in this SO question. The result is the equation of the intersection line in parametric form. Arbitrary values of the parameter t yields a point on this line. Ideally this is what I want to achieve: finding points on the intersection of the planes.

Solving the system of N linear equations with N variables (Ndim) can be done with, for example, Math::Matrix. Given N planes in the form a1x+b1y+c1z+...n1=0 etc., a matrix A is formed with each row being [a1, b1, c1, ..., n1] etc. then: Math::Matrix->new([[a1,b1,c1,...,n1],...])->solve->print;

I have experimented with solve() for the case of M linear equations (Ndim) where M = N-1:

use Math::Matrix;
Math::Matrix->new(
  [1, 2, 1, -1], # ax+by+cz+n=0
  [2, 3, -2, 2]
)->solve->print;
# says
# -7.00000    7.00000 
#  4.00000   -4.00000
[download]

I interpret the result as: the 1st row refers to the 1st dim (aka x) and we have

x -7*t + 7 = 0 => x = 7*t - 7
[download]

The 2nd row yields

y + 4*t -4 = 0 => y = -4*t + 4
[download]

and the missing 3rd dim is the parameter t itself: z = t. And I can get a point on that line by setting t to an arbitrary value. I can verify that this point is on each plane by substituting the point's coordinates into each plane equation.

I have also experimented with transforming the planes matrix to the Reduced Row Echelon Form. Which makes it convenient to work out the "hyperline" equation. This functionality is offered by Math::MatrixLUP :

use Math::MatrixLUP;
my @planes3D = (
  [1, 2, 1, -1], # ax+by+cz+n=0
  [2, 3, -2, 2]
);
print Math::MatrixLUP->new(\@planes3D)->rref
# says
# [1, 0, -7, 7],
# [0, 1, 4, -4]
[download]

Which I interpret as: 1st row has 1 on 1st col, that's the 1st dim (aka x) :

x -7 * t + 7 = 0 => x = 7*t-7
[download]

etc.

So that works well.

Here is a demo I whipped up:

use strict;
use warnings;

use Math::Matrix;
use Math::MatrixLUP;
use Data::Roundtrip qw/perl2dump no-unicode-escape-permanently/;

# this is used for checking if two numbers are sufficiently
# close: abs(x-y)<SMALLNUMBER
# 1E-09 is too small! and for more dimensions 1E-05 is too small
my $SMALLNUMBER = 1E-02;

srand 125; # on the same page

my @planes3D = (
  [1, 2, 1, -1], # ax+by+cz+n=0
  [2, 3, -2, 2]
);

my $sol3D = intersection_by_solve(\@planes3D);
verifysolution(\@planes3D, $sol3D);
$sol3D = intersection_by_rref(\@planes3D);
verifysolution(\@planes3D, $sol3D);

my @planes5D = (
  [1, 2, 1, -2, 3, -1], # ax+by+cz+n=0
  [2, 3, -2, 4, -4, 2],
  [3, -1, -2, -5, 6, 2],
  [-3, 4, 2, 1, 1, 2]
);

my $sol5D = intersection_by_solve(\@planes5D);
verifysolution(\@planes5D, $sol5D);
$sol5D = intersection_by_rref(\@planes5D);
verifysolution(\@planes5D, $sol5D);

# this needs some time ...
#test_in_beast_space();

=pod

# edge cases, 1 dimension is all zeros
# this fails
@planes5D = (
  [0, 2, 1, -2, 3, -1], # ax+by+cz+n=0
  [0, 3, -2, 4, -4, 2],
  [0, -1, -2, -5, 6, 2],
  [0, 4, 2, 1, 1, 2]
);

$sol5D = intersection_by_solve(\@planes5D);
verifysolution(\@planes5D, $sol5D);

$sol5D = intersection_by_rref(\@planes5D);
print Math::MatrixLUP->new(\@planes5D)->rref;
verifysolution(\@planes5D, $sol5D);

# edge cases, 2 dimension is all zeros
# this fails
@planes5D = (
  [0, 0, 1, -2, 3, -1], # ax+by+cz+n=0
  [0, 0, -2, 4, -4, 2],
  [0, 0, -2, -5, 6, 2],
  [0, 0, 2, 1, 1, 2]
);

$sol5D = intersection_by_solve(\@planes5D);
verifysolution(\@planes5D, $sol5D);

$sol5D = intersection_by_rref(\@planes5D);
print Math::MatrixLUP->new(\@planes5D)->rref;
verifysolution(\@planes5D, $sol5D);
=cut

sub intersection_by_solve {
    my $planes = shift;
    die "at least two planes are parallel, can not continue."
       if are_planes_parallel($planes) == 1;
    return Math::Matrix->new($planes)->solve->as_array;
}
sub intersection_by_rref {
    my $planes = shift;
    die "at least two planes are parallel, can not continue."
       if are_planes_parallel($planes) == 1;

    my $sol = Math::MatrixLUP->new($planes)->rref->as_array;
    # leave only the last 2 cols
    return [ map { [ $_->[-2], $_->[-1] ] } @$sol ];
}

sub verifysolution {
    my ($planes, $sol) = @_;
    my $N = scalar(@{$planes->[0]}) - 1;
    #print perl2dump($sol, {terse=>1,indent=>1}). "verifying above sol
+ution in $N dims with ".scalar(@$planes)." planes ...\n";
    print "verifying solution in $N dims with ".scalar(@$planes)." pla
+nes ...\n";
    for (1..1000){
        my $t = -1000 + rand(2000);
        my @point;
        for my $i (0..$#$sol){
            push @point, -($sol->[$i][0] * $t) - $sol->[$i][1];
        }
        push @point, $t; # the last dimension

        # check if point is on all planes
        for my $plane (@$planes){
            my $res = 0.0;
            for my $i (0..$#point){
                $res += $point[$i] * $plane->[$i]
            }
            $res += $plane->[-1]; # the 'n'
            die "point (@point) not in plane (@$plane) ($res!=0) for t
+=$t"
              unless abs($res) < $SMALLNUMBER;
        }
    }
    print "all points lie on the intersecting planes\n";
}

sub verifysolution_by_rref {
    my ($planes, $sol) = @_;
    my $N = scalar(@{$planes->[0]}) - 1;
    for (1..1000){
        my $t = -1000 + rand(2000);
        my @point;
        for my $i (0..$#$sol){
            push @point, -($sol->[$i][0] * $t) - $sol->[$i][1];
        }
        push @point, $t; # the last dimension

        # check if point is on all planes
        for my $plane (@$planes){
            my $res = 0.0;
            for my $i (0..$#point){
                $res += $point[$i] * $plane->[$i]
            }
            $res += $plane->[-1]; # the 'n'
            die "point (@point) not in plane (@$plane) ($res!=0) for t
+=$t"
              unless abs($res) < $SMALLNUMBER;
        }
    }
    print "all points lie on the intersecting planes\n";
}

# Any number of planes, will be checked pairwise
# returns 1 if any two are parallel (and does not check for the rest)
# returns 0 if neither two are parallel (and does not check for the re
+st)
sub are_planes_parallel {
    my $planes = shift;
    my $N = scalar(@{ $planes->[0] }) - 1; # number of dimensions

    # test to see if planes are parallel, in which case there is no in
+tersection:
    # taken pairwise, if the ratio of each of the coefficients (except
+ the constant)
    # is the same then they are parallel
    # special care for zero coefficients, see below:
    for my $i (0..$#$planes){
        my $p1 = $planes->[$i];
        for my $j (($i+1)..$#$planes){
            my $p2 = $planes->[$j];
            my $ratio = undef;
            my $these_two_planes_are_parallel = 1;
            # we should not check the constant (n) in ax+by+cz+n
            INPLANES:
            for my $ii (0..($N-1)){
                # see https://perlmonks.org/?node_id=11160680
                # if both coefficients are zero then it is equivalent
                # to the ratio of these coefficients being the same as
+ the previous ones
                # so we keep checking the other coefficients
                # because that implied planes are "parallel on that di
+mension"

                # two zeroes:
                # they are "parallel to this axis", continue with othe
+r coefficients
                next INPLANES if (abs($p1->[$ii]) < $SMALLNUMBER)
                          && (abs($p2->[$ii]) < $SMALLNUMBER)
                ;

                # one or the other are zero : they are not parallel on
+ this axis
                # and therefore the planes are not parallel,
                # no need to check other coefficients,
                # proceed with other pairs of planes
                if(  (abs($p1->[$ii]) < $SMALLNUMBER)
                 xor (abs($p2->[$ii]) < $SMALLNUMBER) # ^^ is in perl 
+5.40+
                ){
                    # because we complain if is undef below...
                    $ratio = 0 unless defined $ratio;
                    $these_two_planes_are_parallel = 0;
                    last INPLANES
                }                
                unless(defined $ratio){
                    # create a ratio if it was not already created
                    # we will be checking all other ratios to this
                    # if at least one is not the same (wrt range check
+ below)
                    # we break the loop because the planes are not par
+allel
                    $ratio = $p1->[$ii] / $p2->[$ii];
                    next INPLANES;
                }
                my $v = abs($ratio - $p1->[$ii] / $p2->[$ii]);
                if( abs($ratio - ($p1->[$ii] / $p2->[$ii])) > $SMALLNU
+MBER ){
                    $these_two_planes_are_parallel = 0;
                    last INPLANES
                }
            } # INPLANES : for each coefficient of the planes
            # no need to check other planes, at least two are parallel
            if( $these_two_planes_are_parallel > 0 ){
                return 1 # parallel!
            }
        } # end $j plane
    } # end $i plane
    return 0 # not parallel
}

sub test_in_beast_space {
    for (1..1){
        my $N = 666;
        my @planes = (undef) x ($N-1);
        for my $plane (@planes){
            #$plane = [ map { -10 + int rand(20) } (0..$N) ];
            $plane = [ map { -100 + rand(200) } (0..$N) ];
        }
        my $sol = intersection_by_solve(\@planes);
        verifysolution(\@planes, $sol);
        $sol = intersection_by_rref(\@planes);
        verifysolution(\@planes, $sol);
    }
    print "all points lie on the intersecting planes\n";
}
[download]

I have some questions:

Two planes are parallel when their normal vectors (the coefficients of the dimensions: a, b, c, ... (but not the constant n) are multiples. E.g. nv1 = k * nv2. And this translates to checking if all the ratios of the coefficients : a1/a2 = b1/b2 = c1/c2 = ... are the same. My question is: what happens if any coefficient is zero (or actually both coefficients (e.g. a1 and a2) are zero?
How can I calculate the intersection when M < (N-1)? I.e. above I am always checking the intersection of M planes in Ndim where M = N-1. But what if there are even less planes? e.g.
```
my @planes5D = (
  [1, 2, 1, -2, 3, -1], # ax+by+cz+n=0
  [2, 3, -2, 4, -4, 2],
  [3, -1, -2, -5, 6, 2],
  # only 3 planes (it was successful with 4)
)
[download]
```
In short, does the parametric equation of the intersecting "hyperline" contain 2 parameters now?
When I test test_in_beast_space which produces random planes in 666 dimensions, it fails to detect that a point on the intersection "hyperline" lies on all the planes. To do that I substitute the point in each plane equation and expect to have result as zero. However, it's not an exact zero. That's why I am doing a range test to check if it's close to zero. Well the closeness $SMALLNUMBER for 5 dimensions can be as small as 1E-09. But for these many dimensions it can be as low as 1E-02. Is the accuracy lost in summing 666 multiplications that much really?

I guess Math::Matrix::solve() is safe to be used for MxN matrices (e.g. M equations=planes in N unknowns (dimensions) where M<N)? From solve() of Math::Matrix :

Solves a equation system given by the matrix. The number of colums mus
+t be greater than the number of rows. If variables are dependent from
+ each other, the second and all further of the dependent coefficients
+ are 0. This means the method can handle such systems. The method ret
+urns a matrix containing the solutions in its columns or undef in cas
+e of error.
[download]

The edge cases where the 1st dimension is zero for all planes in 5D and when the 1st+2nd are zero, fail. How can I find the intersection in this case?

Edits: updated the demo to correct the sub are_planes_parallel() as per the comments of hippo and LanX about what happens when plane-equation coefficients are zero.

have fun in beast space! bw bliako

Comment on The intersection of M hyperplanes (Ndim) Select or Download Code

Replies are listed 'Best First'.
Re: The intersection of M hyperplanes (Ndim) by hippo (Archbishop) on Jul 18, 2024 at 15:07 UTC
I have also experimented with transforming the planes matrix to the Reduced Row Echelon Form. Isn't it amazing how a simple phrase like "Reduced Row Echelon Form" can transport you back decades to when you last used it (or probably even heard/read it). Thanks for that sudden slip out of the present! My question is: what happens if any coefficient is zero (or actually both coefficients (e.g. a1 and a2) are zero? If both are zero, then that is fine because your normal vectors are not divergent in that axis. If only one or other is zero then they are divergent and the planes are not parallel and will then intersect. You can always rotate your co-ordinate system through something other than π to prove this. That's the only one of your questions I can attempt off the top of my head, sorry. 🦛	[reply]
Re^2: The intersection of M hyperplanes (Ndim) by bliako (Abbot) on Jul 18, 2024 at 15:59 UTC
glad for that flashback So, I understand that if `a1=0 and a2=0` then the vectors are "parallel on that axis", and so I proceed to check the other coefficients in order to determine if the planes are parallel. If one of `a1 OR a2` is zero then I conclude the vectors (and the planes) are not parallel and the test can stop right there.	[reply] [d/l] [select]
Re^3: The intersection of M hyperplanes (Ndim) by LanX (Saint) on Jul 19, 2024 at 09:42 UTC
Look, similar to the rule you already stated: taking two vectors, for them to be collinear there must exist one `k != 0` such that for all coordinate `d: d1 = k * d2` Obviously `d1 == 0 iff d2 == 0` because `k != 0` !!! But if one coordinate is 0 for both it doesn't matter which k is chosen for the other ones. Isn't it obvious now? You check all coordinates in a loop and break if either only one is 0 (next if both are) or the ratio k is not like before (unless undef) Edit FWIW: There is one edge case left, the zero vector. But a normal vector with all coordinates being 0 looks like a bug to me and should raise an exception. Cheers Rolf _{(addicted to the Perl Programming Language :) see Wikisyntax for the Monastery}	[reply] [d/l] [select]
Re^4: The intersection of M hyperplanes (Ndim) by bliako (Abbot) on Jul 19, 2024 at 11:10 UTC
Re^5: The intersection of M hyperplanes (Ndim) by LanX (Saint) on Jul 19, 2024 at 12:58 UTC
Some notes below your chosen depth have not been shown here
Re: The intersection of M hyperplanes (Ndim) by LanX (Saint) on Jul 18, 2024 at 15:11 UTC
> if any coefficient is zero (or actually both coefficients (e.g. a1 and a2) are zero? In order to be collinear `x1==0 <=> x2==0` for all dimensions x Don't fall into the trap to divide by zero. It helps envisioning the 3D case, like `a1 == b1== 0` means the NV1 is on the Z axis. This describes a plane parallel to the XY plane. Cheers Rolf _{(addicted to the Perl Programming Language :) see Wikisyntax for the Monastery}	[reply] [d/l] [select]
Re^2: The intersection of M hyperplanes (Ndim) by bliako (Abbot) on Jul 18, 2024 at 15:23 UTC
`Don't fall into the trap to divide by zero.` hehe computer says no!	[reply] [d/l]
Re: The intersection of M hyperplanes (Ndim) by LanX (Saint) on Jul 26, 2024 at 16:48 UTC
> How can I calculate the intersection when M < (N-1)? I.e. above I am always checking the intersection of M planes in Ndim where M = N-1. But what if there are even less planes? e.g. It depends on your understanding of "calculating", there are multiple ways to represent a subspace. (Update: See Vector space for terminology) The intersection of 3 planes in R_3= R x R x R is obviously a point, the intersection of 2 planes a line. The orthogonal normal vector of a plane is actually a dual space defining the points of the plane by yielding 0 after multiplication (resp a constant, but let's assume all spaces to pass thru the root 0 for simplicity). This goes both ways the plane is the dual space of the line defined by the vector, hence multiplying points of the line with two vectors spawning the plane will yield 0 again. Hence for every subspace of dimension M you can find a dual subspace of dimension K with M+K=N. So in the case of intersecting two planes, their normal vectors can be used to spawn the dual space of the line. That's the plane orthogonal to the line. And it's the trivial condition of intersecting pNV1 = pNV2 = 0, means p is on the line. (I have to complicate all this because you want the abstraction for higher dimensions) Now let's say you use the Gaussian transform for our example. The last line will not yield 1 but 2 values instead of `a . . 0 b . 0 0 c` [download] you get `a . . 0 b c` [download] The line defined by b,c is the projection of the line on the axis plane with x1=0 You can "plot" the intersection by stepping thru x3, to get x2 in the second row and finally get x1 by inserting x2 and x3 in the first row. I hope the abstraction to higher dimensions is clearer now. Dual spaces are relevant because it's easier to make calculations with big M and N, if M-N is small. That's why a simple one dimensional normal vector can define a hyperplane inside a high dimensional space N. Cheers Rolf _{(addicted to the Perl Programming Language :) see Wikisyntax for the Monastery}	[reply] [d/l] [select]
Re^2: The intersection of M hyperplanes (Ndim) by bliako (Abbot) on Jul 30, 2024 at 09:23 UTC
thanks for this detailed explanation. The theory is beyond my level but the last bit (`by stepping thru x3, to get x2 in the second row and finally get x1 by inserting x2 and x3 in the first row`) I understand well. thank you.	[reply] [d/l]
Re^3: The intersection of M hyperplanes (Ndim) by LanX (Saint) on Jul 30, 2024 at 20:44 UTC
Look I know it could be explained better, but without drawings it's quite challenging. 🤷🏻‍♂️ All things involved are handled in the curriculum of multiple school years. Cheers Rolf _{(addicted to the Perl Programming Language :) see Wikisyntax for the Monastery}	[reply]
Re: The intersection of M hyperplanes (Ndim) by LanX (Saint) on Jul 21, 2024 at 00:44 UTC
It would help to know what you're trying to achieve. This looks a lot like a problem from Linear Programming like e.g. the Simplex Algorithm Probably an XY Problem... Cheers Rolf _{(addicted to the Perl Programming Language :) see Wikisyntax for the Monastery}	[reply]
Re^2: The intersection of M hyperplanes (Ndim) by bliako (Abbot) on Jul 24, 2024 at 07:43 UTC
`This looks a lot like a problem from Linear Programming like e.g. the Simplex Algorithm` How? How can the equations I provided be transformed as input to the simplex algorithm (s.a.)? Are there any other objectives other than maximimze or minimize with s.a.? As I don't need either. Also, what happens with s.a. when number of equations is less than Ndim? One day I was contemplating about different supermarkets having different prices for same products and asked "for X money where do you buy more olives and less fetta?"	[reply] [d/l]
Re^3: The intersection of M hyperplanes (Ndim) by LanX (Saint) on Jul 24, 2024 at 12:15 UTC
> "for X money where do you buy more olives and less fetta?" By which criteria ??? Trivially you buy feta in the cheapest supermarket `F` and olives in the cheapest `O` ... Or do you say going to multiple supermarkets isn't practical and you want to minimize the costs for time and fuel? Or do you say you need a minimum amount `f` of feta and `o` of olives and you need the supermarket where you can buy the whole package to the best condition? Cheers Rolf _{(addicted to the Perl Programming Language :) see Wikisyntax for the Monastery}	[reply] [d/l] [select]
Re^4: The intersection of M hyperplanes (Ndim) by bliako (Abbot) on Jul 24, 2024 at 13:00 UTC
Re^5: The intersection of M hyperplanes (Ndim) by LanX (Saint) on Jul 24, 2024 at 17:22 UTC
Some notes below your chosen depth have not been shown here
Re^5: The intersection of M hyperplanes (Ndim) by LanX (Saint) on Jul 24, 2024 at 18:19 UTC