3

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

enter image description here

Community
  • 1
Drey
  • 894
  • 6
  • 11
  • 1
    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$. – Matthew Gunn Mar 10 '17 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 Mar 10 '17 at 19:01
  • 1
    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](http://stats.stackexchange.com/a/34468/919), [using the `contour` function](http://stats.stackexchange.com/questions/24380), and [with a parametric plot](http://stats.stackexchange.com/a/186419/919). – whuber Mar 12 '17 at 22:34

1 Answers1

2

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

Matthew Gunn
  • 20,541
  • 1
  • 47
  • 85
  • Super, thanks. I'm starting to understand. I will try out some things before to accept the answer... – Drey Mar 10 '17 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 Mar 12 '17 at 17:24
  • 1
    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 Mar 12 '17 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 Mar 12 '17 at 18:45
  • 1
    @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). – Matthew Gunn Mar 12 '17 at 18:55
  • Thank you, but do you have an idea how to identify the other points? Seems the rotation alone does not cut it. – Drey Mar 12 '17 at 19:43