# Solving system of linear equations (to determine a boundry)

I'm puzzeled how to programmatically (in R) solve the following linear system:

Given $\mathbf{R} \in \mathbb{R}^{n \times n}$, $\mathbf{R}^{-1}$, and a constant $c$ what is the solution to $\mathbf{u} \in \mathbb{R}^n$ with $\mathbf{u}^T = (u_1, \ldots, u_n)$ for

$\mathbf{u}^T\mathbf{R}^{-1}\mathbf{u} = c$

Lets take the simple case for $n = 2$. Fixing a component, say $u_1$, the solution to $u_2$ can be found by explicitly writing down $u_1(u_1r_{11} + u_2r_{21}) + u_2(u_1r_{12}+u_2r_{22}) - c = 0$ and solving for $u_2$ using quadratic formula. But for more dimensions there must be a better way.

I guess I need to bring it into the form $\mathbf{Ax = b}$, in order to use solve. But I havent figured out yet how exactly.

Right now I'm stuck at the following: let $\mathbf{U} = \left(\begin{matrix} u_1 & \ldots & 0\\ \ldots & \ldots & \ldots\\ 0 & \ldots & 1 \end{matrix}\right)$ and $\mathbf{v}^T = (1, 1, \ldots, u_n)$ with $\mathbf{Uv=u}$ then i would have the fixed terms separated from the variable one ($u_n$) for which I need to determine the value. I can put it into the equation above, but how to proceed? Is this the right way ?

The background is the answer I posted in How to draw confidence areas. I would like to explecitly compute the "exact" threshold boundry. I understand that I need to solve this linear system but I cannot get it quite right yet. I'm unsatisfied with the two possible solutions: 1. using Quadratic formula to hard code the solution and 2. using optimize routine. The first one would only work for 2 dimensions and the second one would be unreliable (because upto two different solutions are possible for every x).

Furthermore, I think there should be a concise solution.

edit (12.03) Thank you for the response. I played with the solution but still have some question.

So as far as I understood, compute_scale would compute my decision boundry. Since I have two possibilities for $\gamma$, i.e. positive and negative, I can compute the critical values. However, if I plot them, I only get the half truth. I tinkered, but havent figured out how to compute the complete boundry. Any advice?

compute_stat <- function(v, rmat) {
transv <- qnorm(v)
return(as.numeric(transv %*% rmat %*% transv))
}

compute_scale <- function(v, rmat) {
gammavar <- sqrt(threshold / (v %*% rmat %*% v))
return(c(pos = pnorm(v * gammavar), neg = pnorm(v * (-gammavar))))
}

Rg <- matrix(c(1, .1, .2, 1), ncol = 2)#matrix(c(1,.01,.99,1), ncol = 2)
Rginv <- MASS::ginv(Rg)

gridval <- seq(10^-2, 1 - 10^-2, length.out = 100)
thedata <- expand.grid(x = gridval,
y = gridval)

thestat <- apply(thedata, 1, compute_stat, rmat = Rginv)

threshold <- qchisq(1 - 0.8, df = 2)
colors <- ifelse(thestat < threshold, "#FF000077", "#00FF0013")
#png("boundry2.png", 640, 480)
plot(y ~ x, data = thedata, bg = colors, pch = 21, col = "#00000000")

theboundry <- t(apply(thedata, 1, compute_scale, rmat = Rginv))
points(pos1 ~ pos2, data = theboundry, col = "blue")
points(neg1 ~ neg2, data = theboundry, col = "purple")
#dev.off()

• You don't have a system of equations. You have a single equation. Let $A = R^{-1}$. Your single equation is $\sum_{i=1}^n \sum_{j=1}^n u_i u_j A_{ij} = c$. Commented Mar 10, 2017 at 17:03
• Please specify what form a solution ought to take, given that generically an $n-1$ dimensional manifold of solutions exists. There are usually two ways to draw pictures of such solutions: implicitly, by giving formulas they satisfy (which is what you already have), and parametrically, by giving $n-1$ functions of one real variable that collectively trace out the solutions. But others are available, such as a method to produce a new solution from any given one (which thereby can be iterated). A "solution" can also be conceptual, too: in your application it is simply a rescaled hypersphere.
– whuber
Commented Mar 10, 2017 at 19:01
• Since you haven't responded to my requests for clarification, about all I can do is refer to you to various worked examples and let you pick among them: using the ellipse function, using the contour function, and with a parametric plot.
– whuber
Commented Mar 12, 2017 at 22:34

I understand your problem to be given an $n$ by $n$ matrix $R$ and scalar $c$, find a vector $\mathbf{u}$ such that $\mathbf{u}'R^{-1}\mathbf{u}=c$.

First observe:

• You have $n$ unknowns (since $\mathbf{u}$ is an $n$ by 1 vector)
• $\mathbf{u}'R^{-1}\mathbf{u}=c$ is a single equation. (It isn't a system of equations.)

In general, there won't be a unique solution $\mathbf{u}$. Almost any vector will work if it is properly scaled.

### Solution:

Pick some arbitrary vector $\mathbf{a}$. Let $\mathbf{u} = \lambda \mathbf{a}$. Then $\mathbf{u}'R^{-1}\mathbf{u}=c$ becomes $\lambda^2 \mathbf{a}'R^{-1}\mathbf{a} = c$. Solving for the scalar $\lambda$ we have $\lambda = \sqrt{\frac{c}{\mathbf{a}'R^{-1}\mathbf{a}}}$.

For any vector $\mathbf{a}$ such that $\mathbf{a}'R^{-1}\mathbf{a} \neq 0$, we'll have the solution:

$$\mathbf{u} = \lambda \mathbf{a}\quad \text{where} \quad \lambda = \sqrt{\frac{c}{\mathbf{a}'R^{-1}\mathbf{a}}}$$

• Super, thanks. I'm starting to understand. I will try out some things before to accept the answer...
– Drey
Commented Mar 10, 2017 at 17:34
• I added some insights, but I haven't figured out how to compute the complete boundry. Would you have any advice? Thanks!
– Drey
Commented Mar 12, 2017 at 17:24
• The final formula for $\lambda$ has to be understood as a generalized square root: both the positive and negative roots need to be included. Otherwise you will get only half of the solutions.
– whuber
Commented Mar 12, 2017 at 17:32
• thats is exactly what I did @whuber. I included the postive and negative roots in the code. Actually one value of $\lambda$ only covers quarter of the solution.
– Drey
Commented Mar 12, 2017 at 18:45
• @Drey You don't need the negative lambda even if you iterate over vectors $\mathbf{a}$ that point in all directions. From your graph, it appears you only have vectors with an angle from 0 to 90 degrees, so that gets you the boundary that's in the first quadrant (and the negation gets you the boundary in the 3rd quadrant). Commented Mar 12, 2017 at 18:55