finding the orthogonal vectors for a series of vectors

Hello,

I'm looking for some help to calculate the orthogonal vectors to a series of vectors. I have vector xhat which is (about700000x3) and every line or row of it is the vector that I want to calculate the two orthogonals for. I am aware of the null function however I'm struggling to implement it for the data type type I have. Any help is appriciated.

[Torsten]

I am aware of the null function however I'm struggling to implement it for the data type type I have.

The null function is the easiest way to achieve what you want. I don't understand your comment about the data type.

a = [1 3 -6];
w = null(a)

w = 3×2

-0.4423 0.8847 0.8295 0.3410 0.3410 0.3180

a*w 

ans = 1×2

1.0e-14 * -0.1221 0.1776

Answer 1

[Bjorn Gustavsson]

This is how I would do it:

nVecs = 12
allVecs = randn(nVecs,3);
for i1 = size(allVecs,1):-1:1
  u_test = randn(1,3);
  while is_parallel(allVecs(i1,:),u_test) % You'll have to write this test yourself
    u_test = randn(1,3);
  end
  % now we have a vector u_test that are not parallel to allVec(i1,:)
  % First perpendicular vector:
  VecPerp1(i1,:) = cross(allVecs(i1,:),u_test);
  % then a second perpendicular to both allVecs(i1,:) and VecPerp1(i1,:) is
  % the cross-product of those two.
  VecPerp2(i1,:) = cross(allVecs(i1,:),VecPerp1(i1,:));
end

You will have to clean up this snippet a bit with normalizzations etc. Perhaps you also want some other characteristics of the perpendicular vectors - but you can always make some pair of linear combinations or rotations of them afterwards, or select the test-vector u_test a bit differently. This fullfills the conditions you presented.

HTH

[Pouya]

Thank you Bjorn,

I understand the majority of the code you suggested. One question that I have regarding your suggested answer is the u_test vector. Can you please let me know what it is and how I should create it?

Edit: I combined your answer with Dyuman's below and made it work. Thank you both!

[Bjorn Gustavsson]

You're welcome, happy that we helped.

My algorithm is based of the idea that a random array, the u_test, is not likely to be parallel to the current vector in allVecs. If it is we just pick a new one. Then we cross-multiply that arbitrary vector (which also are unlikely to be perpendicular to the currenct vector) with the current vector. That gives us a vector that should be perpendicular to the current vector. Then we can create a vector that is perpendicular to both the current vector in allVecs and that first perpendicular vector with another cross-multiplication.

[Bruno Luong]

I see many of you want to generate random vector and make sure it is not parallel to the original vector, say u.

One non-rejection way doing that is

r = zeros(1,3);
[~,i] = min(abs(u));
r(i) = 1;

such r is never parallel to u.

EDIT: I just saw John d'Errico proposes the same idea.

[John D'Errico]

@Bruno Luong - surely it is that great minds think alike.

[Bjorn Gustavsson]

 @Bruno Luong, well I was lazy and thought that a random vector in 3-D would be sufficiently non-parallel to get away with this, then I tossed in a check just to be safe, because of said lazyness I couldn't think up your general step (in my work I have alway had a safe non-parallel direction in [0 0 1] for this task). I've checked how likely one would be to have to reject the first random vector - and overnight my laptop did some 1e9-1e10 tests without a rejection, so in general one might get away with that, but it would be an absollutely terrible programming habit.

[David Goodmanson]

Hi Bjorn,

What you are calling an absolutely terrible programming habit, I don't know that I agree. Let a target unit vector and a randomly chosen unit vector differ by some amount (the difference being a vector that is almost exactly perpendicular to each of them). Suppose that if the difference is >= than, say, 1e-7 then the vectors are distinguishable. That's a reasonable tolerance. With uniform spacing on average, there are 4e14 such vectors on the unit sphere that exceed the tolerance. One could wait for a very long time before a randomly selected vector didn't work.

[Bjorn Gustavsson]

@David Goodmanson, Yeah, maybe that was putting it a bit too bluntly. Since Bruno's and your solutions aren't that demanding they are preferable - and I have had a habit of sailing a bit too close to the edge and gotten bitten a couple of times. But I'm happy to get upgraded to "not ideal programming habit" or however far you allow it.

Answer 2

[Dyuman Joshi]

Another method will be -

for i=1:size(vec,1)
    v=vec(i,:);
    %generate a random vector
    r=rand(1,3);
    %check it is not parallel to original vector
    while ~cross(r,v)
        r=rand(1,3);
    end
    %first perpendicular vector
    q1=r-dot(r,v)/dot(v,v)*v;
    %second perpendicular vector
    q2=cross(q1,v);
end
%store the vectors as you like

There are some trivial methods as well -

%for example let your vector be [a b c]
%then perpendicular vectors are - 
q1 = [-b a 0];
q2 = [-c 0 a];
q3 = [0 c -b];

[Pouya]

Hi Dyuman,

I think your first suggestion is working but the only problem is q1 and q2 that I'm getting are 1x3 when they are supposed to be size(vec)x3. Any suggestions?

[Dyuman Joshi]

So, as I mentioned in a comment in my code, it's upto you how you want to store them.

You can pre-allocate empty/zeros erray then use indexes or concatenate them. Like -

%pre-allocation
q1=zeros(size(vec));
q2=zeros(size(vec));
for i=1:size(vec,1)
    ... %skipping the code
    p=r-dot(r,v)/dot(v,v)*v;
    q1(i,:)=p;
    q2(i,:)=cross(q1(i,:),v);
end

Answer 3

[John D'Errico]

I think you are confused. Really, this is not a difficult problem. You need to understand the linear algebra of it, or you will still be confused. You have an array of size 700000x3, or whatever. For each row of this array, you want to find TWO vectors normal to the corresponding row. Of course, they will not be unique, since we can always rotate them arbitrarily around the axis of the original vector.

Lacking any data from you, I'll create some completely random vectors.

V = randn(700000,3); % My personal set of initial vectors

We will want that set to have unit norm though, as that makes the later computations easier.

n = size(V,1); % how many vectors are there?
Vnormalized = normalize(V,2,'norm'); % scales every row to have unit 2-norm 

Next, for each vector, we need to find a new vector that will NEVER be parallel to the vectors in Vnormalized. I'll be tricky here, to insure we will always succeed. The idea will be to find the element with the smallest absolute value in each vector, and modify that element. That will create a new vector that is ABSOLUTELY assured to be not parallel.

[~,ind] = min(abs(Vnormalized),[],2);
Q1 = eye(3);
Q1 = normalize(Q1(ind,:) + Vnormalized,2,'norm'); % normalize Q1

But while Q1 is not parallel to V, it is not orthogonal either. We can make that happen now using a dot product, subtracting off the component of Q1 that is still parallel to Vnormalized.

Q1 = Q1 - sum(Q1.*Vnormalized,2).*Vnormalized; % requires R2016b or later
Q1 = normalize(Q1,2,'norm'); % re-normalize Q1

Ok. We now have TWO vectors, Vnormalized and Q1, bioth of which have unit norm, and which are orthogonal. The third is easy. Just use a cross product. And since I've already set up Vnormalized and Q1 to be unit vectors, then so will be the cros product.

Q2 = cross(Vnormalized,Q1);

The resulting rows of Q1 and Q2 are now each unit norm, and they are orthogonal to the corresponding rows of V. For example, consider the first row of V.

M = [Vnormalized(1,:);Q1(1,:);Q2(1,:)]

M = 3×3

0.4339 -0.0649 -0.8986 0.0282 0.9979 -0.0585 0.9005 -0.0000 0.4348

This next result will be the identity matrix if we were successful, to within floating point trash in the least significant bits.

format long g
M'*M

ans = 3×3

1 1.36670842986521e-16 0 1.36670842986521e-16 1 -2.91736562051799e-16 0 -2.91736562051799e-16 1

And that is as good as we can do.

Answer 4

[Bruno Luong]

If you have recent MATLAB (from 2021b)

% Invent some test data
xhat = rand(10,3);

[B,~,~]=pagesvd(reshape(xhat',3,1,[]));
N=B(:,2:3,:); 
% N = pagetranspose(N) if prefer 2 row vectors for each page

% Check, each of two columns of N(:,:,k) is orthogonal to xhat(k,:) for all k
xhat(1,:)*N(:,:,1) % almost 0

ans = 1×2

1.0e-16 * -0.2776 -0.5551

xhat(end,:)*N(:,:,end) % almost 0

ans = 1×2

1.0e-15 * 0.0009 0.1110

[Pouya]

Bruno, I can't thank you enough! This worked perfectly!

[Bruno Luong]

you are welcome

Answer 5

[Bruno Luong]

A direct method for vectors in R^3

a=rand(100,3); % but be non zeros

% Generate b and c orthogonal to a
a2 = a(:,[2 3 1]); a3 = a(:,[3 1 2]);
b = a2-a3;
c= (a2+a3)-(sum(a,2).^2./sum(a.^2,2)-1).*a;
% Normalize a, b, c if required ...

% Check orthogonality
norm(sum(a.*b,2),Inf)

ans = 1.1102e-16

norm(sum(a.*c,2),Inf)

ans = 1.6098e-15

norm(sum(c.*b,2),Inf)

ans = 1.6653e-16

There is some spectial treatment needs to be applied when a has 3 component identical (parallel to [1,1,1]). But I'll leave this detail out, unless someone requrires it.

[Bruno Luong]

"There is some spectial treatment needs to be applied when a has 3 component identical"

What the heck I post it here:

degenerate = all(b==0,2);
p = sum(degenerate);
b(degenerate,:) = repmat([0 1 -1], [p,1]);
c(degenerate,:) = repmat([-2 1 1], [p,1]);

Answer 6

[David Goodmanson]

Hi Pouya,

Here is another method, somewhat related to Bruno's

a = 2*rand(1000,3)-1;

m1 = tril(ones(3,3)) -triu(ones(3,3))
m2 = [1 1 1]'*[-1 1 -1];
m3 = [1 -1 1]'*[1 1 1];
b = a*m1;
c = (a.*a)*m2 + a.*(a*m3); 

% check
max(abs(diag(a*b')))
max(abs(diag(b*c')))
max(abs(diag(c*a')))

[Bruno Luong]

And just like me, your method needs to pay attention when a is parallel to [1 -1 1].

[David Goodmanson]

Yes. I imagine that any method with a linear relationship between a and b is going to have a vector in 'a' that doesn't work.

<added>

If a is nx1 and M is nxn, then if

b = M*a  and  a'*b = 0   for all nonzero a,

then M must be antisymmetric. In an odd number of dimensions, M cannot be full rank, so null(M) is the required vector in 'a' such that M*a = 0.

[Bruno Luong]

can you prove there is a non trivial a such that c=0, assuming only

b = a*m1;
c = (a.*a)*m2 + a.*(a*m3); 
b'*a = 0
c'*a = 0

and not a specific m1, m2, m3 you provide?

[David Goodmanson]

Hi Bruno,

do you mean b and c as defined, with m1, m2, m3 such that

a~=0, b~=0, c~=0,  a*b'=0,  b*c'=0,  c*a'=0

almost always, then to find an exceptional a such that

b~=0, c=0,

I don't think that is possible. Are you relaxing the requirement that

b*c'=0 ? 

[Bruno Luong]

Yes I forgot b'*c = 0.

What I mean is finding a such that : a~=0 and c=0.

The question is simple, you have shown that there is non trivial a that makes b=m1*a=0. I ask the same questiin with c instead of b.

If you tackle such prove (or disprove), you can use b (or not) and the fact that b'*c = 0, or a'*m1*c = 0, meaning a and c are conjugate via m1 an anti-symmetric matrix, and singular for odd dimension. I'm not sure it helps.

Obviously in your specific example a = (1,-1,0) works for both.

[David Goodmanson]

Hi Bruno,

the expressions I came up with are basically

b = a x (1,-1,1)      and
c =  a x b  = a ((1,-1,1).a) - (1,-1,1) (a.a) 

rigged up to operate on an nx3 matrix of 3d row vectors. So it isn't possible for c to equal 0 unless b also equals 0.

I have spent too much time on the following problem. Given a vector a = (x,y,z), find another vector linearly independent of it (after that one can get an orthogonal set using cross products). It's easy enough to do if you allow an if-then-else procedure to take care of at least one special case, but I was trying to do so in such a way that all special cases, such as (1,0,0), (x,y,0) , (+-1,+-1,+-1) etc. are handled automatically. A fixed point theorem may come into play, so I am allowing noncontinuous functions such as abs, mod, rem. Square roots with specified sign, trig functions, whatever. I'm still not there.

[Bruno Luong]

@David Goodmanson As we discuss earlier with John for a given not-zero u vector

v = zeros(size(u));
[~,i] = min(abs(u));
v(i) = 1;

v is always independent to u.

If one prefer an MATLAB one-line expression

v = accumarray(find(abs(u)==min(abs(u)),1),1,size(u))

NOTE: Too bad matlab doesn't have build in argmin function on discrete data (return a second output of min on vector).

The expression can be shorter and more readable

v = accumarray(argmin(abs(u)),1,size(u))

[David Goodmanson]

Thanks Bruno. Yes, how could I have been so blind?