For a physics application, I need to find the solution $X \in \mathbb R^{n \times m}$ to the matrix equation
$$ AX=B $$
where $A \in \mathbb R^{n \times n}$ and $B \in \mathbb R^{n \times m}$, and $n \gg m$. Roughly speaking, $n = 2 \cdot 10^4$ and $m = 150$.
Partitioning the matrices into blocks,
$$ \begin{pmatrix}I_m & A_1 \\ 0 & A_2 \end{pmatrix} \begin{pmatrix}X_1 \\ X_2 \end{pmatrix} = \begin{pmatrix}B_1 \\ B_2 \end{pmatrix}, $$
where $A_2$, $X_1$, and $B_1$ are square. I am only interested in obtaining $X_1$, but it is necessary to solve for the entire $X$ in order to have a fully determined system. $X_1$ should be symmetric, and I am looking for a way to solve this equation in such a way the forces the symmetry of $X_1$.
Currently, I use an LAPACK routine which solves the system for each column of $X$, and the resulting $X_1$ is symmetric up to some precision, but I am having issues with numerical instability so I am hoping that knowing $X_1$ should be symmetric might add some extra condition to the system that improves the stability of the solution.
Is there a way to rewrite this expression in such a way that when I solve it $X_1$ is guaranteed to be exactly symmetric?