Finding a vector orthogonal to n n-dimensional vectors

Question

First off, I would like to state that I was wondering about where I should have posted this question. Math SE or StackOverflow? But it seemed that this question was more of a math question than a programming one; if it's not the case, please pardon me.

I'm trying to find a vector (not all vectors; only one would suffice) perpendicular to n other n-dimensional vectors. e.g. I have 3 three-dimensional vectors and I'd like to find an extra one which is orthogonal to all the other 3. I'm trying to implement this program in C#, using Cramer's Rule for matrix calculation. My algorithm is based on the notion that the dot product of the desired orthogonal vector with all the other vectors should be 0: Basically, if our vectors are $v_1, v_2, .., v_n$ and our desired vector is $k$ then $dot(k,v_1) = 0, dot(k,v_2)=0,...,dot(k,v_n)=0$. In my program, I generate 3 three-dimensional vectors with random component numbers and try to find vector $k$ orthogonal to all 3 of them. The chances are I'm gonna face something like this:

$v_1= (1,3,6)$, $v_2 = (4, 9, -1)$, $v_3 = (-3,1,12)$

These vectors at the top are the randomly-generated ones. Then the program I'm developing, plugs-in the data into a system of linear equations' matrix: $$ \begin{bmatrix} 1 & 3 & 6 \\ 4 & 9 & -1 \\ -3 & 1 & 12 \\ \end{bmatrix}\begin{bmatrix} x \\ y \\ z\\ \end{bmatrix}=\begin{bmatrix} 0 \\ 0\\ 0\\ \end{bmatrix} $$ and I will get $x=y=z=0$ which basically means $k=(x,y,z)=(0,0,0)$ which is expected, because there's no non-zero vector orthogonal to all three of them. It's all decent, until the random vector generators produce $v_1=(1,2,3),v_2=(2,4,6),v_3=(0,-1,2)$. Which would give me: $$ \begin{bmatrix} 1 & 2 & 3 \\ 2 & 4 & 6 \\ 0 & -1 & 2 \\ \end{bmatrix}\begin{bmatrix} x \\ y \\ z\\ \end{bmatrix}=\begin{bmatrix} 0 \\ 0\\ 0\\ \end{bmatrix} $$ There actually are an infinite number of vectors orthogonal to all three of $v_1,v_2$ and $v_3$ but with my current algorithm, I can't calculate even one of them, because the determinant for the coefficient matrix is zero. I can plug in $x=1$ and find $y$ and $z$, but that only works handsomely for 3-dimensional vectors. I want to develop an algorithm that works for any n-dimensional vector. Hence, n could be 10, and therefore I can't go around plugging in random values for $x,y,z,w,..$ because that would not work. I think the approach I've taken for this program is very likely wrong. I took a look at the Gram-Schmidt Process, but I'm not sure if that will help me. All I want is calculating one, nice, cute orthogonal vector out of all vectors orthogonal to n n-dimensional vectors. Any help would be greatly appreciated. I'm adding C# code (in UnityEngine) for what I'm trying to achieve, in case it helps.

CramersRule.cs:

using System;
using System.Collections.Generic;
using static System.Linq.Enumerable;
using UnityEngine;

public static class CramersRule
{

public static int[] SolveCramer(int[][] equations)
{
    int size = equations.Length;
    if (equations.Any(eq => eq.Length != size + 1)) throw new ArgumentException($"Each equation must have {size + 1} terms.");
    int[,] matrix = new int[size, size];
    int[] column = new int[size];
    for (int r = 0; r < size; r++)
    {
        column[r] = equations[r][size];
        for (int c = 0; c < size; c++)
        {
            matrix[r, c] = equations[r][c];
        }
    }
    return Solve(new SubMatrix(matrix, column));
}

private static int[] Solve(SubMatrix matrix)
{
    int det = matrix.Det();
    if (det == 0)
    {
        return new int[] { };
    }

    int[] answer = new int[matrix.Size];
    for (int i = 0; i < matrix.Size; i++)
    {
        matrix.ColumnIndex = i;
        answer[i] = matrix.Det() / det;
    }
    return answer;
}

//Extension method from library.
static string DelimitWith<T>(this IEnumerable<T> source, string separator = " ") =>
    string.Join(separator ?? " ", source ?? Empty<T>());

private class SubMatrix
{
    private int[,] source;
    private SubMatrix prev;
    private int[] replaceColumn;

    public SubMatrix(int[,] source, int[] replaceColumn)
    {
        this.source = source;
        this.replaceColumn = replaceColumn;
        this.prev = null;
        this.ColumnIndex = -1;
        Size = replaceColumn.Length;
    }

    private SubMatrix(SubMatrix prev, int deletedColumnIndex = -1)
    {
        this.source = null;
        this.prev = prev;
        this.ColumnIndex = deletedColumnIndex;
        Size = prev.Size - 1;
    }

    public int ColumnIndex { get; set; }
    public int Size { get; }

    public int this[int row, int column]
    {
        get
        {
            if (source != null) return column == ColumnIndex ? replaceColumn[row] : source[row, column];
            return prev[row + 1, column < ColumnIndex ? column : column + 1];
        }
    }

    public int Det()
    {
        if (Size == 1) return this[0, 0];
        if (Size == 2) return this[0, 0] * this[1, 1] - this[0, 1] * this[1, 0];
        SubMatrix m = new SubMatrix(this);
        int det = 0;
        int sign = 1;
        for (int c = 0; c < Size; c++)
        {
            m.ColumnIndex = c;
            int d = m.Det();
            det += this[0, c] * d * sign;
            sign = -sign;
        }
        return det;
    }

    public void Print()
    {
        for (int r = 0; r < Size; r++)
        {
            Debug.Log(Range(0, Size).Select(c => this[r, c]).DelimitWith(", "));
        }
    }
}

}

SceneManager.cs:

using System.Collections;
using System.Collections.Generic;
using UnityEngine;

public class SceneManager : MonoBehaviour
{

void Start()
{
    //Code gets called on initialization, same as main(void) in C or main(String[] args) in Java
    if (false)
    {
        InitRandomVectors();
    } else
    {

        //Vector 0
        List<int> vector0 = new List<int>(); vector0.Add(3); vector0.Add(-1); vector0.Add(4);
        List<int> vector1 = new List<int>(); vector1.Add(-2); vector1.Add(3); vector1.Add(1);
        List<int> vector2 = new List<int>(); vector2.Add(-4); vector2.Add(6); vector2.Add(2);

        vectors.Add(vector0); vectors.Add(vector1); vectors.Add(vector2); //This is gonna result in a 0 determinant

    }
    InitCramersRulesCoeffResultMatrix();
    DisplayRandomValsAndSolution();
}

private void DisplayRandomValsAndSolution()
{
    int[] solution = CramersRule.SolveCramer(matrices);
    if (solution.Length == 0)
    {
        //Faced a 0 det
        return;
    }
    for (int i = 0; i < vectorCount; i++)
    {
        int[] vector = new int[vectorDimension];
        for (int j = 0; j < vector.Length; j++)
        {
            vector[j] = vectors[i][j];
        }
        Debug.Log($"Vector {i} = " + DisplayAsVector(vector));
    }
    Debug.Log($"Solution Vector = {DisplayAsVector(solution)}");
}

private string DisplayAsVector(int[] vectorComponents)
{
    string s = "(";
    int vectorDimension = vectorComponents.Length;
    for (int j = 0; j < vectorDimension; j++)
    {
        s += ((j != 0) ? ", " : "") + vectorComponents[j];
    }
    s += ")";
    return s;
}

private int[][] matrices = null;

public void InitCramersRulesCoeffResultMatrix()
{
    matrices = new int[vectorCount][];
    for (int i = 0; i < vectorCount; i++)
    {
        int[] vectorCoefficientsWith0Result = new int[vectorDimension + 1];
        for (int j = 0; j < vectorDimension; j++)
        {
            vectorCoefficientsWith0Result[j] = vectors[i][j];
        }
        vectorCoefficientsWith0Result[vectorDimension] = 0;
        matrices[i] = vectorCoefficientsWith0Result;
    }
}

private const int vectorDimension = 3;
private const int vectorCount = 3;
private const int vectorComponentAbsoluteRange = 5;

private List<List<int>> vectors = new List<List<int>>();

public void InitRandomVectors()
{
    for (int i = 0; i < vectorCount; i++)
    {
        List<int> innerVectorComponents = new List<int>();
        for (int j = 0; j < vectorDimension; j++)
        {
            int randomComponentValue = Random.Range(-vectorComponentAbsoluteRange, +vectorComponentAbsoluteRange);
            innerVectorComponents.Add(randomComponentValue);
        }
        vectors.Add(innerVectorComponents);
    }
}

}

I thought maybe I would use diophantine equations (with my very limited knowledge on the subject), but that would make it one big hassle to implement for any n-dimensional vector. So I passed on that.

Answer 1

all you need to do is to do ordinary row reduction, where the requirement be that you have correctly implemented rational numbers. Only elementary operations with rational numbers will be allowed. In the end, you will have a rational vector that is a null vector of the original matrix, or a reduced matrix, upper triangular with all nonzero entries on the diagonal, in which case it cannot be done.

Answer 2

If you put the matrix in reduced row echelon form, you can put anything in to all the components that don’t have a leading one in their column of the matrix, then each component with a leading one in its column gets the negative dot product of the rest of the row the leading one is from, with the things you put in to the empty rows.

Answer 3

You are trying to find an elements of $S^\perp$ with $S = \langle v_1, \dots, v_n\rangle$. By linear algebra (assuming we are working in $\mathbb{R}^n$) we know that

$$ \dim S^\perp = \dim \mathbb{R}-\dim S = n-\dim S. $$

Thus, if $v_1,\dots,v_n$ are linearly independent, they form a basis for $S$ and thus $\dim S = n, \dim S^\perp = 0$. This shows that the only possible choice for such a vector is $0$.

Now, suppose that $v_1,\dots,v_n$ are not linearly independent. In particular, the matrix

$$ A = \begin{pmatrix} \mid & \mid & \cdots & \mid \\ v_1 & v_2 & \cdots & v_n\\ \mid & \mid & \cdots & \mid \end{pmatrix} $$

cannot be invertible, meaning that we have some nonzero solution for

$$ Ax = 0, $$

and therefore there exists a vector $x \in \mathbb{R}^n$ such that $\langle v_i,x\rangle = 0$ for all $1 \leq i \leq n$ as you desire.

Cramer's rule will only work when $A$ is invertible, and in that case you already know that the only possible solution for you problem is to take the zero vector.

When $A$ is not invertible, the vectors you are looking for are the solutions to $Ax = 0$. This can be computed via Gauß-Jordan elimination, for example (i.e. taking $A$ to is row-echelon form). Most likely $C\#$ has packages which do this already.