# finding the orthogonal vectors for a series of vectors

111 views (last 30 days)
Pouya on 22 Jul 2022
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 on 22 Jul 2022
Edited: Torsten on 22 Jul 2022
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

Bjorn Gustavsson on 22 Jul 2022
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
Bjorn Gustavsson on 5 Aug 2022
@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.

Dyuman Joshi on 22 Jul 2022
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];
##### 2 CommentsShowHide 1 older comment
Dyuman Joshi on 23 Jul 2022
So, as I mentioned in a comment in my code, it's upto you how you want to store them.
You can initialise/pre-allocate empty/zeros erray then use indexes or concatenate them. Like -
%concatenation
q1=[];
q2=[];
for i=1:size(vec,1)
... %skipping the code
p=r-dot(r,v)/dot(v,v)*v;
q1=[q1;p];
q2=[q2;cross(p,v)];
end
%indexing
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

John D'Errico on 22 Jul 2022
Edited: John D'Errico on 22 Jul 2022
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.

Bruno Luong on 22 Jul 2022
Edited: Bruno Luong on 22 Jul 2022
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
##### 2 CommentsShowHide 1 older comment
Bruno Luong on 25 Jul 2022
you are welcome

Bruno Luong on 23 Jul 2022
Edited: Bruno Luong on 23 Jul 2022
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 on 23 Jul 2022
"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]);

David Goodmanson on 25 Jul 2022
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')))
David Goodmanson on 4 Aug 2022
Thanks Bruno. Yes, how could I have been so blind?