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
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
The 2nd row yields
y + 4*t -4 = 0 => y = -4*t + 4
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]
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
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";
}
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)
)
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.
- 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