How to generate a randomly oriented, high dimensio

Question

I want to generate a randomly oriented circle in R^n. I have been able to successfully generate points on the surface of an n-sphere. I read that you can intersect it with a plane and get a circle, but I don't know how to do this in Python.

Or is there any other way to generate it in Python?

Thanks!

Answer 1

As a complete python implementation :

N=3 # dimensions
R=1 # radius
npoints=1000
from numpy.random import rand
from numpy.linalg import norm
from numpy import dot,sin,cos,outer,pi
V1,V2,origin=rand(3,N)
u1=V1/norm(V1)
u2=V2/norm(V2)
V3=u1-dot(u1,u2)*u2
u3=V3/norm(V3)
print(norm(u2),norm(u3),dot(u2,u3))
#1.0 1.0 0.0
#u2,u3 is orthonormed 
theta=np.arange(0,2*pi,2*pi/npoints)

circle=origin+R*(outer(cos(theta),u2)+outer(sin(theta),u3))

import matplotlib as mpl
from mpl_toolkits.mplot3d import Axes3D
import matplotlib.pyplot as plt
fig = plt.figure()
ax = fig.gca(projection='3d')
ax.plot(*circle.T)
plt.show()

Answer 2

You need to use 2 basis vectors for this.

1. Generate Random unit N-D vector U

   so just create array of N random numbers ui = <-1,+1> and then normalize to unit size by dividing them all by the size sqrt(u0^2+u1^2+...) Beware The U must not be zero !!!

2. Generate N-D vector V which is perpendicular to U

   In 2D and 3D is this simple as you can just swap x,y and negate one in 2D or use cross product in 3D. Unfortunately The cross product is not defined in N-D (at least not to my knowledge) So you need to generate non zero vector for which dot product is zero:

   (U.V) = dot(U,V) = 0
   u0.v0 + u1.v1 . u2.v2 + ... = 0
   

   So you can create some elements as random and fit the rest (so the sum is zero). To equalize you can use abs higher vi where ui is abs lower... After this normalize V to unit size.

   In the C++ example bellow I used this approach:

   1. compute random v[i]=< -1 , +1 > such that partial dot product is lovering in magnitude (so just compute sign so u[i]*v[i] is negative to partial dot product)
   2. find any non zero u[i] and recompute v[i] so dot product is zero
   3. normalize V to unit size

   Beware also V must not be zero !!!

3. Circle in N-D

   use simple parametric equation:

   P(t) = center + radius*cos(t)*U + radius*sin(t)*V
   

   Where t = < 0 , 2*Pi >

Here simple C++ example of generating U,V:

//---------------------------------------------------------------------------
const int n=5; // dimension count
//---------------------------------------------------------------------------
double nd_dot(double *a,double *b)  // = (a.b)
    {
    int i; double c;
    for (c=0.0,i=0;i<n;i++) c+=a[i]*b[i];
    return c;
    }
//---------------------------------------------------------------------------
void nd_unit(double *a) // a = a / |a|
    {
    int i; double c;
    for (c=0.0,i=0;i<n;i++) c+=a[i]*a[i];
    c=sqrt(c); if (c>1e-10) c=1.0/c; else c=0.0;
    for (i=0;i<n;i++) a[i]*=c;
    }
//---------------------------------------------------------------------------
double nd_length(double *a) // = |a|
    {
    int i; double c;
    for (c=0.0,i=0;i<n;i++) c+=a[i]*a[i];
    return sqrt(c);
    }
//---------------------------------------------------------------------------
void nd_rnd(double *a)  // a = (? ? ? ... ?) , |a| = 1.0
    {
    int i; double c;
    for (;;)
        {
        for (i=0;i<n;i++) a[i]=Random()-0.5;
        for (c=0.0,i=0;i<n;i++) c+=a[i]*a[i];
        if (c>1e-20) break;
        }
    c=1.0/sqrt(c);
    for (i=0;i<n;i++) a[i]*=c;
    }
//---------------------------------------------------------------------------
void nd_perpendicular(double *a,double *b)  // a = (? ? ? ... ?) , |a| = 1.0 , (a.b) = 0.0
    {
    if (n==1)   // 1D case
        {
        if (fabs(b[0]>1e-10)) a[0]=1.0; else a[0]=1.0;
        return;
        }
    int i; double c1,s;
    for (c1=0.0,i=0;i<n;i++)                    // c1 = dot(a,b)
        {
        s=1.0;                                  // s = signum of a[i] so (a.b) -> 0
        if (b[i]<0.0) s=-1.0;
        if (c1*s>0.0) s=-s;
        a[i]=s*Random();
        c1+=a[i]*b[i];
        }
    for (i=0;i<n;i++)                           // correct single element so (a.b) = 0
     if (fabs(b[i])>1e-10)
        {
        c1-=a[i]*b[i];
        a[i]=-c1/b[i];
        break;
        }
    nd_unit(a);
    }
//---------------------------------------------------------------------------
AnsiString nd_prn(double *a) // this is for printing you can ignore this
    {
    int i; AnsiString s="( ";
    for (i=0;i<n;i++) s+=AnsiString().sprintf("%6.3lf ",a[i]);
    s+=" )"; return s;
    }
//---------------------------------------------------------------------------
// This is usage example in VCL you can ignore this:
__fastcall TForm1::TForm1(TComponent* Owner):TForm(Owner)
    {
    double u[n],v[n];
    Randomize();
    for (int i=0;i<10;i++)
        {
        nd_rnd(u);
        nd_perpendicular(v,u);
        mm_log->Lines->Add("U = "+nd_prn(u));
        mm_log->Lines->Add("V = "+nd_prn(v));
        mm_log->Lines->Add("(U.V) = "+AnsiString().sprintf("%6.3lf",nd_dot(u,v)));
        }
    }
//-------------------------------------------------------------------------

And the result of this for confirmation of approach:

U = ( -0.526 -0.623 -0.535 -0.215 -0.044  )
V = (  0.620 -0.023 -0.287 -0.682 -0.260  )
(U.V) = -0.000
U = (  0.444  0.282  0.517 -0.203  0.644  )
V = (  0.072 -0.564 -0.499  0.134  0.640  )
(U.V) = -0.000
U = (  0.636  0.339  0.357 -0.022 -0.594  )
V = ( -0.559 -0.015 -0.108 -0.497 -0.655  )
(U.V) =  0.000
U = ( -0.626 -0.195  0.491 -0.374 -0.435  )
V = ( -0.722 -0.063 -0.541  0.039  0.424  )
(U.V) =  0.000
U = ( -0.532 -0.481  0.467 -0.517  0.019  )
V = (  0.536 -0.716 -0.344 -0.205 -0.199  )
(U.V) = -0.000
U = (  0.696 -0.588  0.018 -0.078  0.405  )
V = ( -0.106 -0.514 -0.645 -0.079 -0.550  )
(U.V) =  0.000
U = (  0.072  0.679  0.124 -0.204 -0.690  )
V = (  0.995 -0.058  0.041  0.060  0.037  )
(U.V) = -0.000
U = ( -0.742 -0.283  0.579 -0.150  0.109  )
V = ( -0.043 -0.798 -0.512 -0.308 -0.066  )
(U.V) = -0.000
U = (  0.606  0.389 -0.551  0.041 -0.420  )
V = ( -0.457 -0.153 -0.489 -0.691 -0.228  )
(U.V) =  0.000
U = (  0.947 -0.225  0.156 -0.075  0.151  )
V = (  0.125 -0.153 -0.043 -0.034 -0.979  )
(U.V) = -0.000

[Edit1] randomness issue

As @RoobieNuby pointed out method above can have randomness issues (the planes will not have uniform distribution) After some testing I confirmed it:

So I tested 3 and 2 basis vectors form and the difference can be visually seen. On the right is the above method in 5D (truncated to 2D view) and on the left is new form obtained like this:

//$$---- Form CPP ----
//---------------------------------------------------------------------------
#include <vcl.h>
#include <math.h>
#pragma hdrstop
#include "Unit1.h"
//---------------------------------------------------------------------------
#pragma package(smart_init)
#pragma resource "*.dfm"
TForm1 *Form1;
Graphics::TBitmap *bmp=NULL;
int xs=0,ys=0,**pyx=NULL;

const int n=5;
//---------------------------------------------------------------------------
double nd_dot(double *a,double *b)  // = (a.b)
    {
    int i; double c;
    for (c=0.0,i=0;i<n;i++) c+=a[i]*b[i];
    return c;
    }
//---------------------------------------------------------------------------
void nd_unit(double *a) // a = a / |a|
    {
    int i; double c;
    for (c=0.0,i=0;i<n;i++) c+=a[i]*a[i];
    c=sqrt(c); if (c>1e-10) c=1.0/c; else c=0.0;
    for (i=0;i<n;i++) a[i]*=c;
    }
//---------------------------------------------------------------------------
double nd_length(double *a) // = |a|
    {
    int i; double c;
    for (c=0.0,i=0;i<n;i++) c+=a[i]*a[i];
    return sqrt(c);
    }
//---------------------------------------------------------------------------
void nd_rnd(double *a)  // a = (? ? ? ... ?) , |a| = 1.0
    {
    int i; double c;
    for (;;)
        {
        for (i=0;i<n;i++) a[i]=Random()-0.5;
        for (c=0.0,i=0;i<n;i++) c+=a[i]*a[i];
        if (c>1e-20) break;
        }
    c=1.0/sqrt(c);
    for (i=0;i<n;i++) a[i]*=c;
    }
//---------------------------------------------------------------------------
void nd_perpendicular(double *a,double *b)  // a = (? ? ? ... ?) , |a| = 1.0 , (a.b) = 0.0
    {
    if (n==1)   // 1D case
        {
        if (fabs(b[0]>1e-10)) a[0]=1.0; else a[0]=1.0;
        return;
        }
    int i; double c,s;
    for (c=0.0,i=0;i<n;i++)                     // c = dot(a,b)
        {
        s=1.0;                                  // s = signum of a[i] so (a.b) -> 0
        if (b[i]<0.0) s=-1.0;
        if (c*s>0.0) s=-s;
        a[i]=s*Random();
        c+=a[i]*b[i];
        }
    for (i=0;i<n;i++)                           // correct single element so (a.b) = 0
     if (fabs(b[i])>1e-10)
        {
        c-=a[i]*b[i];
        a[i]=-c/b[i];
        break;
        }
     nd_unit(a);
    }
//---------------------------------------------------------------------------
void nd_perpendicular(double *a,double *b,double *c) // a = (? ? ? ... ?) , |a| = 1.0 , (a.b) = 0.0 , (a.c) = 0.0
    { // this is not in-place so (&a != &b) and (&a != &c)
    int i,e; double ab,ac;
    // spec cases
    if (n<3) { for (i=0;i<n;i++) a[i]=0.0; return; }
    // init
    for (i=0;i<n;i++) a[i]=1.0;
    ab = nd_dot(a,b);
    ac = nd_dot(a,c);
    // tune until dot products near zero
    for (e=0;(e<1000)&&(fabs(ab)+fabs(ac)>=1e-5);e++)   // max 1000 iterations so it will not hang up
     for (i=0;i<n;i++)
        {
        ab-=a[i]*b[i];
        ac-=a[i]*c[i];
             if (fabs(b[i])>fabs(c[i])) a[i]=-ab/b[i];
        else if (fabs(b[i])<fabs(c[i])) a[i]=-ac/c[i];
        else if (fabs(ab)>=fabs(ac)) a[i]=-ab/b[i];
        else                         a[i]=-ac/c[i];
        ab+=a[i]*b[i];
        ac+=a[i]*c[i];
        }
    nd_unit(a);
    }
//---------------------------------------------------------------------------
AnsiString nd_prn(double *a)
    {
    int i; AnsiString s="( ";
    for (i=0;i<n;i++) s+=AnsiString().sprintf("%6.3lf ",a[i]);
    s+=" )"; return s;
    }
//---------------------------------------------------------------------------
// VCL application init (can ignore this)
__fastcall TForm1::TForm1(TComponent* Owner):TForm(Owner)
    {
    bmp=new Graphics::TBitmap;
    bmp->HandleType=bmDIB;
    bmp->PixelFormat=pf32bit;
    if (bmp==NULL) Application->Terminate();
    }
//-------------------------------------------------------------------------
// VCL application resize event (can ignore this)
void __fastcall TForm1::FormResize(TObject *Sender)
    {
    if (bmp==NULL) return;
    xs=ClientWidth-mm_log->Width;
    ys=ClientHeight;
    bmp->SetSize(xs,ys);
    xs=bmp->Width;
    ys=bmp->Height;
    if (pyx) delete[] pyx;
    pyx=new int*[ys];
    for (int y=0;y<ys;y++) pyx[y]=(int*)bmp->ScanLine[y];
    Paint();
    }
//---------------------------------------------------------------------------
// VCL application exit (can ignore this)
void __fastcall TForm1::FormDestroy(TObject *Sender)
    {
    if (bmp) delete   bmp; bmp=NULL;
    if (pyx) delete[] pyx; pyx=NULL;
    }
//---------------------------------------------------------------------------
// VCL application repaint (consider this void main()...)
void __fastcall TForm1::FormPaint(TObject *Sender)
    {
    if (bmp==NULL) return;
    if (pyx==NULL) return;
    int i,e,x,y;
    double u[n],v[n],w[n];
    double a,da,x0,y0,r,c,s;
    Randomize();

    // clear screen (xs,ys is resolution, pyx is direct pixel access to bitmap bmp)
    for (y=0;y<ys;y++)
     for (x=0;x<xs;x++)
      pyx[y][x]=0x00000000;

    da=1.0*M_PI/180.0;              // angle step
    x0=0.5*xs;                      // center
    y0=0.5*ys;
    r=150.0;                        // radius
    // 100 random circles
    for (i=0;i<100;i++)
        {
        // 3 vector form
        nd_rnd(w);                  // W random unit vector (normal of the plane)
        nd_perpendicular(v,w);      // V perpendicular to W
        nd_perpendicular(u,v,w);    // U perpendicular to V and W
        // old 2 vector form
//      nd_rnd(u);
//      nd_perpendicular(v,u);

        for (e=1,a=0.0;e;a+=da)     // circle points loop
            {
            if (a>=2.0*M_PI) { e=0; a=0.0; }
            c=r*cos(a);
            s=r*sin(a);
            x=double(x0+(c*u[0])+(s*v[0]));     // use just x,y for render
            y=double(y0+(c*u[1])+(s*v[1]));
            if ((x>=0)&&(x<xs)&&(y>=0)&&(y<ys)) // if inside screen
             pyx[y][x]=0x00FFFFFF;              // render it
            }
        // debug info log to see the U,V,W as numbers
        if (i==0)
            {
            mm_log->Text="";
            mm_log->Lines->Add("U = "+nd_prn(u));
            mm_log->Lines->Add("V = "+nd_prn(v));
            mm_log->Lines->Add("W = "+nd_prn(w));
            mm_log->Lines->Add("(U.V) = "+AnsiString().sprintf("%6.3lf",nd_dot(u,v)));
            mm_log->Lines->Add("(U.W) = "+AnsiString().sprintf("%6.3lf",nd_dot(u,w)));
            mm_log->Lines->Add("(W.V) = "+AnsiString().sprintf("%6.3lf",nd_dot(w,v)));
            }
        }
    // refresh Application window with bitmap (bacbuffering)
    Canvas->Draw(0,0,bmp);
    }
//---------------------------------------------------------------------------
// VCLo mouse wheel event ... just force repaint on mouse wheel to check the randomness in time (can ignore this)
void __fastcall TForm1::FormMouseWheel(TObject *Sender, TShiftState Shift, int WheelDelta, TPoint &MousePos, bool &Handled)
    {
    Paint();
    }
//---------------------------------------------------------------------------

The important stuff here is void nd_perpendicular(double *a,double *b,double *c); function which returns perpendicular vector to both b,c. To make this robust (In case randomness matters) in N-D You should have n vectors not just 3.

The operation of function is similar to the previous one. The difference is only in that we optimize more that one dot product at once. So you need to choose which dot product to optimize per axis (based on abs coordinate value of the perpendicular vector). So you should rewrite this to optimize n-1 dot products at once.

1. create 1st vector (normal)
2. create 2nd vector (1 dot product to optimize)
3. create 3th vector (2 dot products to optimize)
4. create 4th vector (3 dot products to optimize)
5. ...
6. create n-th vector (n-1 dot products to optimize)
7. use any 2 (except the first) vectors as basis (you can select them randomly)
8. produce your circle

It looks like just 3 vectors are enough for n>=3 (so you can ignore #4,#5) but to confirm that for sure some sort of statistical analysis would be needed which I am too lazy to do.

[Edit2] vector perpendicular to more vectors ...

Have implemented the dot product perpendicular vector computation from the oher answer. And computation of perpendicular vector to multiple vectors:

//---------------------------------------------------------------------------
void nd_perpendicular(double *a,double *b)  // a = (? ? ? ... ?) , |a| = 1.0 , (a.b) = 0.0
    {
    int i; double tmp[n],dot,len;
    // |dot| limit to consider as non parallel
    len=0.7*nd_length(b);
    // tmp = random unit vector not parallel to b
    for (;;)
        {
        nd_rnd(tmp);
        dot=nd_dot(b,tmp);
        if (fabs(dot)<=len) break;
        }
    // make a unit and perpendicular to b
    for (i=0;i<n;i++) a[i]=tmp[i]-(dot*b[i]);
    nd_unit(a);
    }
//---------------------------------------------------------------------------
void nd_perpendicular(double *a,double *v0,double *v1,double *v2=NULL,double *v3=NULL,double *v4=NULL,double *v5=NULL,double *v6=NULL,double *v7=NULL,double *v8=NULL,double *v9=NULL)
    {
    // this is not in-place so (&a != &b) and (&a != &c)
    // a = unit vector perpendicular to all not NULL vectors v0,v1,v2,v3,...
    const int M=10; // vi operands max count
    int i,j,k,e,m; double dot[M],*v[M]={v0,v1,v2,v3,v4,v5,v6,v7,v8,v9},q;
    // spec cases
    if (n<3) { for (i=0;i<n;i++) a[i]=0.0; return; }
    // init
    for (i=0;i<n;i++) a[i]=1.0;
    nd_rnd(a);
    for (m=0,j=0;j<M;j++) if (v[j]) { m=j+1; dot[j]=nd_dot(a,v[j]); } else dot[j]=0.0;
    // tune until dot products near zero
    for (e=0;e<1000;e++)    // max 1000 iterations so it will not hang up
        {
        for (q=0.0,j=0;j<m;j++) q+=fabs(dot[j]);
        if (q<1e-3) { e=-1; break; }
        // k = index of abs max dot
        for (k=0,j=0;j<m;j++) if (fabs(dot[j])>=fabs(dot[k])) k=j;
        // i = index of abs max coordinate of v[k]
        for (i=0,j=0;j<n;j++) if (fabs(v[k][j])>=fabs(v[k][i])) i=j;
        // update dots and a[i]
        for (j=0;j<m;j++) dot[j]-=a[i]*v[j][i];
        a[i]=-dot[k]/v[k][i];
        for (j=0;j<m;j++) dot[j]+=a[i]*v[j][i];
        }

    if (e>=0)
     for (e=0;e<1000;e++)   // max 1000 iterations so it will not hang up
        {
        for (q=0.0,j=0;j<m;j++) q+=fabs(dot[j]);
        if (q<1e-3) { e=-1; break; }
        for (i=0;i<n;i++)
            {
            // k = index of abs max vec
            for (k=0,j=0;j<m;j++) if (fabs(v[j][i])>=fabs(v[k][i])) k=j;
            // update dots and a[i]
            for (j=0;j<m;j++) dot[j]-=a[i]*v[j][i];
            a[i]=-dot[k]/v[k][i];
            for (j=0;j<m;j++) dot[j]+=a[i]*v[j][i];
            }
        }
    nd_unit(a);
    }
//---------------------------------------------------------------------------

Too boost robustness It uses 2 different iterative approaches. Here example of usage (n=5):

// nd_perpendicular test
double V[n][n]; int x,y;
nd_rnd(V[0]);
nd_perpendicular(V[1],V[0]);
nd_perpendicular(V[2],V[1],V[0]);
nd_perpendicular(V[3],V[2],V[1],V[0]);
nd_perpendicular(V[4],V[3],V[2],V[1],V[0]);
for (x=0;x<n;x++)
 for (y=x+1;y<n;y++)
  mm_log->Lines->Add(AnsiString().sprintf("(V%i.V%i) = %6.3lf",x,y,nd_dot(V[x],V[y])));

Answer 3

The following is assuming that by high dimensional circle you mean the intersection of an N-sphere with an N-hyperplane containing the sphere's centre (i.e. an N-1-sphere in N+1-space).

You can do the following (assuming you are in N+1-space, considering the N sphere):

• pick a normal vector n, your circle will be the subset of the sphere that is orthogonal to n
• use SVD to get an orthonormal basis B for the orthogonal complement of n
• use your method to create points p on the N-1 sphere and multiply them by B (B is N+1 x N, so Bp will be in N+1-space)

Answer 4

Since you already know how to generate a random point on the surface of an n-sphere, just generate two such points, call them P1 and P2. These will determine the plane in which the circle will lie.

I am assuming that both these points are a distance of 1 from the origin, (which will be true if you picked 1 as the radius of your n-sphere). If not, divide each point by its length.

Now we want to make P2 perpendicular to P1. This can be done by

P2 = P2 - dot(P1, P2) P1
P2 = P2 / || P2 ||

Then for each point, you can generate a random angle between 0 and 2pi. Convert this angle to a point on the circle by:

cos(angle) * P1 + sin(angle) * P2.