How to formulate intersection point of n-D equations with each having distinct origin? Let say we have $n$-dimensional line/plane/hyperplane $n$ represented by $n$ amount of equations. Note that each of the equations is crossing different origin $O\in\mathbb{R}^n$.
What would be the simplest way to:

*

*compute, and optionally

*verify (type, existence) of their intersection $I\in\mathbb{R}^n$?

A simplest example perhaps would be like this:
$$\begin{equation}\begin{aligned}
x_1 + y_1 &= a_1\\
\end{aligned}\end{equation}\tag{1}\label{eq1}$$
$$\begin{equation}\begin{aligned}
x_2 + y_2 &= a_2\\
\end{aligned}\end{equation}\tag{2}\label{eq2}$$
with \eqref{eq1} is having origin $O_1$ and \eqref{eq2} is having origin of $O_2$.
Plotted like this:

2 lines with different origin in $\mathbb{R}^2$
Let say $I$ is graphically proven to exists like above.
What is the analytical solution to reach $I_{x}$ and $I_{y}$?
Plus point if it can be formulated as numerically stable algebraic problem.
 A: If we use vector notation, point $\vec{p}$ is on hyperplane $\vec{n}, d$ if and only if
$$\vec{p} \cdot \vec{n} = d \tag{1}\label{1}$$
But, if the hyperplane definition is with respect to origin $\vec{o}$, then it is
$$(\vec{p} - \vec{o}) \cdot \vec{n} = d \tag{2a}\label{2a}$$
which is equivalent to
$$\vec{p} \cdot \vec{n} = d + \vec{o} \cdot \vec{n} \tag{2b}\label{2b}$$
In other words, having the hyperplane be defined with respect to some point not at origin, only changes the scalar component (signed distance) of the hyperplane.  In absolute coordinates, the scalar component is increased by the dot product between the hyperplane normal and the point used as the origin for the hyperplane definition.
If we use $\vec{p} = (p_1 , \dots , p_N)$ for the coordinates of a point to be considered, $\vec{n} = (n_1 , \dots , p_N)$ for the hyperplane normal, and $\vec{o} = (o_1, \dots, o_N)$ for the point used as the origin for the definition of the hyperplane, we can write $\eqref{2a}$ as
$$(p_1 - o_1) n_1 + \dots + (p_N - o_N) n_N = d \tag{3a}\label{3a}$$
and $\eqref{2b}$ as
$$p_1 n_1 + \dots + p_N n_N = d + o_1 n_1 + \dots + o_N n_N \tag{3b}\label{3b}$$
or, in sum form, as
$$\sum_{i=1}^{N} p_i n_i = d + \sum_{i=1}^N o_i n_i \tag{3c}\label{3c}$$
Let's say you have two hyperplanes, $(n_1, \dots, n_N ; d)$ and $(u_1, \dots, u_N ; p)$, with "origins" $(o_1, \dots, o_N)$ and $(g_1, \dots, g_N)$:
$$\left\lbrace ~ \begin{aligned}
\sum_{i=1}^N p_i n_i &= d + \sum_{i=1}^N o_i n_i \\
\sum_{i=1}^N p_i u_i &= p + \sum_{i=1}^N g_i u_i \\
\end{aligned} \right . \tag{4} \label{4}$$
where both equations in $\eqref{4}$ are only true for $(p_1, \dots, p_N)$ at the intersection of the two hyperplanes.
The only difference to the standard situation is the extra constant sum added to the scalar component, depending only on the "origin" (with respect to which the hyperplane was defined) and on the hyperplane normal.

We have $N$ hyperplanes.  Let's use notation $\vec{n}_i = (n_{i, 1}, \dots, n_{i, N})$ for their normal vectors, and $d_i$ for their signed distances from absolute origin, i.e. including the constant sum on the right side in $\eqref{3c}$ or $\eqref{4}$.
If we form a matrix $\mathbf{M}$, with each row consisting of a hyperplane normal,
$$\mathbf{M} = \left[ \begin{matrix}
n_{1, 1} & n_{1, 2} & \dots & n_{1, N-1} & n_{1, N} \\
n_{2, 1} & n_{2, 2} & \dots & n_{2, N-1} & n_{2, N} \\
\vdots & \vdots & ~ & \vdots & \vdots \\
n_{N-1, 1} & n_{N-1, 2} & \dots & n_{N-1, N-1} & n_{N-1, N} \\
n_{N, 1} & n_{N, 2} & \dots & n_{N, N-1} & n_{N, N} \\
\end{matrix} \right ] \tag{5}\label{5}$$
and a column vector $y$ from the signed distances including the constant "origin" compensation sums,
$$y = \left[ \begin{matrix}
d_1 \\
d_2 \\
\vdots \\
d_{N-1} \\
d_{N} \\
\end{matrix} \right ] \tag{6}\label{6}$$
with $x$ being the column vector for the intersection point, the system of equations becomes
$$\mathbf{M} x = y \tag{7a}\label{7a}$$
This has a solution if $\mathbf{M}$ is invertible,
$$x = \mathbf{M}^{-1} y \tag{7b}\label{7b}$$
In a computer program, $\mathbf{M}$ can be near-singular if some of the normals are almost linearly dependent (i.e, a pairwise dot product is close to the product of their Euclidean norms; or, equivalently, the normals are almost parallel or almost opposite).  Obviously, if two or more of the normals are parallel or opposite, then the intersection is no longer a point.
Unless you have unstated additional requirements, basically any linear algebra library should work for this just fine.
