Nearest intersection point to many lines c#

Question

Based on the good work done [here][1] I am converting a simplified version of the python code to C#. Below is what I have so far (Sorry - I am new to c#). I am stuck on the last statement least squares equivalent in c#. Unfortunately it appears NumSharp hasn't incorporated least square code yet. The inputs to my least squares statement is a R matrix and q vector. I tried using the MathNet library but I can't get it to work with the inputs I have. Is there a better way in C# to compute this least squares? Thanks.

Ref: nearest intersection point to many lines in python

using MathNet.Numerics.LinearAlgebra;
using MathNet.Numerics.LinearRegression;
using NumSharp;
/*
//Python code below:
import numpy as np

def LS_intersect(p0,a0):
    """
    :param p0 : Nx3 (x,y) position coordinates, z=1
    :param a0 : angles in degrees for each point in p0
    :return: least squares intersection point of N lines from eq. 13 in 
             http://cal.cs.illinois.edu/~johannes/research/LS_line_intersect.pdf
    """    

    # create direction vectors with magnitude = 1
    n = []
    for a in a0:
        n.append([np.cos(np.radians(a)), np.sin(np.radians(a))])

    n = np.array(n) #cos & sin of vectors
    
    pos = p0[:,0:2] #vector (x,y) IP coordinates
    
    # generate the array of all projectors 
    nnT = np.array([np.outer(nn,nn) for nn in n ]) 
    ImnnT = np.eye(len(pos[0]))-nnT # orthocomplement projectors to n

    # now generate R matrix and q vector
    R = np.sum(ImnnT,axis=0)
    q = np.sum(np.array([np.dot(m,x) for m,x in zip(ImnnT,pos)]),axis=0)

    # and solve the least squares problem for the intersection point p 
    return np.linalg.lstsq(R,q,rcond=None)[0]

#sample data 
pa = np.array([[-10.0,  10.0, 1. ],
       [-18.0,  -7.0,  1. ],
       [-15.0,  -15.0,  1. ]])
paa = [-45.0,  0.0, 45.0]

solution = LS_intersect(pa,paa)
print(solution)
#Answer should be: [ 5.44009282e-15 -3.50000000e+00]

Start C# code:
*/
// See https://aka.ms/new-console-template for more information
Console.WriteLine("Hello, World!");
//import numpy as np

//double LS_intersect(float a0)

//sample data 
float[,] p0 = new float[,] { { -10.0f, 10.0f }, { -18.0f, -7.0f }, { -15.0f, -15.0f } };
float[] a0 = new float[] { -35.0f, 5.0f, 40.0f };

double[,] v = new double[a0.Length, 2];
for (int i = 0; i < a0.Length; i++)
{
    v[i, 0] = Convert.ToDouble(Math.Cos(ConvertToRadians(a0[i])));
    v[i, 1] = Convert.ToDouble(Math.Sin(ConvertToRadians(a0[i])));
}
double[,] n = (double[,])np.array(v); //array of cos & sin of direction vectors

// generate the array of all projectors

//nnT = np.array(np.outer(nn,nn) for nn in n );

double[][,] nnT = new double[n.GetLength(0)][,];

double[,] nn = new double[1, 2];
for (int i = 0; i < n.GetLength(0); i++)
{
    nn[0, 0] = n[i, 0];
    nn[0, 1] = n[i, 1];
    nnT[i] = (double[,])np.array(np.outer(nn, nn));
}

// ImnnT = np.eye(len(pos[0]))-nnT # orthocomplement projectors to n
double[][,] ImnnT = new double[n.GetLength(0)][,];
for (int i = 0; i < n.GetLength(0); i++)
{
    ImnnT[i] = (double[,])np.subtract(np.eye(2), np.array(nnT[i])); // orthocomplement projectors to n
}

// now generate R matrix and q vector
//    R = np.sum(ImnnT,axis=0)

double[,] R = new double[2, 2];
for (int i = 0; i < n.GetLength(0); i++)
{
    R = (double[,])np.add(R, ImnnT[i]);
}

//var q = np.sum(np.array(np.dot(m,x) for m,x in z.zip(ImnnT,p0)),0);

double[,] qp = new double[n.GetLength(0), 2];
double[] p0e = new double[2];
for (int i = 0; i < n.GetLength(0); i++)
{
    p0e[0] = p0[i, 0];
    p0e[1] = p0[i, 1];
    Vector<double> V = CreateVector.DenseOfArray(p0e);
    Matrix<double> M = CreateMatrix.DenseOfArray(ImnnT[i]);

    Vector<double> MV = M * V;

    qp[i, 0] = MV[0];
    qp[i, 1] = MV[1];
}

double[] q = new double[2];
for (int i = 0; i < n.GetLength(0); i++)
{
    q[0] = q[0] + qp[i, 0];
    q[1] = q[1] + qp[i, 1];
}

double ConvertToRadians(double angle)
{
    return (Math.PI / 180) * angle;
}

//double[] answer = new double[2];
double[][] RR = new double[2][];
for (int i = 0; i < 2; i++)
{
    double[] row = new double[2];
    for (int j = 0; j < 2; j++)
    {
        row[j] = R[i, j];
    }
    RR[i] = row;
}

double[] p = MultipleRegression.QR(   // ********problem is in the statement
    new[] { new[] { 1.0, 4.0 }, new[] { 2.0, 5.0 } },
    new[] { 15.0, 20 },
   intercept: true);

//var answer = MultipleRegression.QR(RR, q, intercept: true);

//var solution = LS_intersect(p0, a0);

// answer should be: [ 3.8085419  -1.94293867]