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")

enter image description here

  • 1
    $\begingroup$ 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$. $\endgroup$ Commented Mar 10, 2017 at 17:03
  • $\begingroup$ 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. $\endgroup$
    – whuber
    Commented Mar 10, 2017 at 19:01
  • 1
    $\begingroup$ 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. $\endgroup$
    – whuber
    Commented Mar 12, 2017 at 22:34

1 Answer 1


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.


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}}}$$

  • $\begingroup$ Super, thanks. I'm starting to understand. I will try out some things before to accept the answer... $\endgroup$
    – Drey
    Commented Mar 10, 2017 at 17:34
  • $\begingroup$ I added some insights, but I haven't figured out how to compute the complete boundry. Would you have any advice? Thanks! $\endgroup$
    – Drey
    Commented Mar 12, 2017 at 17:24
  • 1
    $\begingroup$ 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. $\endgroup$
    – whuber
    Commented Mar 12, 2017 at 17:32
  • $\begingroup$ 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. $\endgroup$
    – Drey
    Commented Mar 12, 2017 at 18:45
  • 1
    $\begingroup$ @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). $\endgroup$ Commented Mar 12, 2017 at 18:55

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.