User brian borchers - MathOverflow

Answer by Brian Borchers for solve non-convex quadratic constrained quadratic programming

2013-05-21T03:21:15Z

It's helpful if you cite the paper in which you saw something that you're asking a question about- we could provide a better answer if we knew where the question came from.

First, assume without loss of generality that $A$ and $C$ are symmetric matrices. It's easy to take these quadratic forms and write them in terms of symmetric matrices.

I believe that this problem transformation requires that $A$ be a positive semidefinite matrix- see below.

We'll begin with your second problem and show that it is equivalent to the original problem. We begin with

$\min \mbox{tr}(AB) $

subject to

$\mbox{tr}(CB)=1$

$ \beta \geq 0 $

$\left[ \begin{array}{cc} 1 & \beta^{T} \\ \beta & B \end{array} \right] \succeq 0 $

By Schur's theorem, the constraint

$\left[ \begin{array}{cc} 1 & \beta^{T} \\ \beta & B \end{array} \right] \succeq 0 $

is equivalent to

$B \succeq \beta \beta^{T} $

Note that we've implicitly restricted $B$ to being a symmetric matrix.

Next, write $B$ as

$B= \beta \beta^{T} + LL^{T} $

where $L$ is the (slightly generalized) Cholesky factor of $B- \beta \beta^{T}$. If $B-\beta \beta^{T}$ is singular, then $L$ would be singular or even $0$.

Then our problem is equivalent to

$\min \mbox{tr}( A ( \beta \beta^{T} + LL^{T} ) ) $

subject to

$\mbox{tr}(CB)=1$

$ \beta \geq 0 $

$ B=\beta \beta^{T} + LL^{T} $

By the property $\mbox{tr}(DEF)=\mbox{tr}(FDE)$, this problem is equivalent to

$\min \mbox{tr}(\beta^{T}A\beta + L^{T}AL) $

subject to

$\mbox{tr}(CB)=1$

$ \beta \geq 0 $

$ B=\beta \beta^{T} + LL^{T} $

Note that because $A \succeq 0$, $\mbox{tr}(L^{T}AL) \geq 0$ for all $L$, and the term is minimized when $L=0$.

Thus our problem is equivalent to:

$\min \mbox{tr}(\beta^{T}A\beta) $

subject to

$\mbox{tr}(CB)=1$

$ \beta \geq 0 $

$ B=\beta \beta^{T} $

Since $B=\beta \beta^{T}$, and by the cyclic property of $\mbox{tr}()$, this is equivalent to

$\min \mbox{tr}(\beta^{T}A\beta) $

subject to

$\beta^{T}C\beta=1$

$ \beta \geq 0 $

This was your original problem.

Answer by Brian Borchers for What does "Vertex Solution" mean?

2013-05-20T13:47:48Z

In the context of linear programming, and assuming that you're using the simplex method to solve your LP's rather than an interior point method, it's most likely that the author means "basic feasible solution" (BFS) here. In geometrical terms, the basic feasible solutions of an LP are vertices of the polytope of feasible solutions.

The number of vertices of a polytope defined by a finite system of linear equalitions and inequalities is finite and bounded by a function involving the number of variables and constraints.

Since each iteration of the simplex method ends with a basic feasible solution, any optimal solution returned by the simplex method will be a BFS. An algorithm which solves a sequence of linear programming problems in which the constraints do not change (but the objective function varies) by using the simplex method will produce a sequence of basic feasible solutions from this finite set.

Answer by Brian Borchers for minimization of a function when the feasible set is an unbounded cone

2013-03-24T04:16:35Z

No. The problem with this simple minded approach is that the sequence of constraints that you add might go on forever without adding a critical constraint.

Consider the following example problem.

$\min x$

subject to

$ x \ge 0$

$ x \ge -1-1/n, \;\;\; n=1, 2, ...$

Now, suppose that you start with the inequality $x \ge -2$, and that in each successive iteration you find that the next inequality of the form $x \ge -1-1/n$ is violated and add it to the problem. Your sequence of solutions will have $x$ converging to -1, which clearly isn't optimal.

You'll need to specify additional structure in the problem and/or use a more sophisticated algorithm.

Answer by Brian Borchers for Projection and Positive matrices

2012-12-26T21:06:03Z

What do you mean by "positive linear combination"?

Try G=I. If you select any column of G, you'll find that its projection onto the space spanned by the other columns of G is the 0 vector.

Answer by Brian Borchers for Relating the angle between two vectors to max and min eigenvalues

2012-12-26T03:32:37Z

Here's a restatement of the problem for those who don't want to find the paper referenced in the question.

We're interested in a symmetric and positive definitive matrix $\Sigma$ that has been orthogonally diagonalized as

$ \Sigma=Q^{T}\mbox{diag}(\theta_{1}^{2},...,\theta_{n}^{2})Q $

where $\theta_{1}^{2} \geq \theta_{2}^{2} \geq ... \geq \theta_{n}^{2} > 0$, and we let $\theta_{\max}^{2}=\theta_{1}^{2}$ and $\theta_{\min}^{2}=\theta_{n}^{2}$.

Also, let

$\Theta^{2}=\mbox{diag}(\theta_{1}^{2},...,\theta_{n}^{2})$.

We want to show for any nonzero vector $\alpha$,

$ \frac{\alpha^{T}\Sigma^{-1}\alpha}{\sqrt{\alpha^{T}\alpha}\sqrt{\alpha^{T}\Sigma^{-2}\alpha}} \geq \frac{\theta_{\max}\theta_{\min}}{(\theta_{\max}^{2}+\theta_{\min}^{2})/2} $

Note that you can assume without loss of generality that $\alpha$ is of length 1. (just scale the length of alpha out of everything on the left hand side of the inequaility.)

Also, by using the substitution $x=Q\alpha$, you can reduce this to a problem about the digonal matrix $\Theta$, and then reduce the matrix-vector products to sums. We then want to show for all vectors $x$ of length 1,

$ \frac{\sum_{i=1}^{n} x_{i}^{2}\theta_{i}^{-2}} {\sqrt{\sum_{i=1}^{n}x_{i}^{2}\theta_{i}^{-4}}} \geq \frac{\theta_{\max}\theta_{\min}}{(\theta_{\max}^{2}+\theta_{\min}^{2})/2} $

Unfortunately, I don't see any easy way to proceed from here.

Answer by Brian Borchers for How to solve a system of linear equations without storing the matrix?

2012-12-07T14:21:30Z

Being able to get elements of the matrix isn't very useful (particularly if you don't know where the nonzero elements of the matrix are without checking.)

Iterative methods can be useful if you have the ability compute matrix vector products $Mx$. You haven't said whether this is possible.

It seems quite likely that there's some special structure to your particular problem that would make it possible to simplify this computation. You haven't told us anything about where the system of equations comes from- perhaps if you explained this in some detail we could suggest ways to proceed.

Answer by Brian Borchers for What software one needs to solve a big linear system on a small computer?

2012-10-19T13:40:02Z

Linear systems of 10,000 equations in 10,000 unknowns can easily be solved in a few seconds using double precision floating point arithmetic on typical consumer grade PC's and even laptop computers.

If you're writing your main program in a compiled language like C or Fortran, you consider using the LAPACK and BLAS libraries to do this kind of work.

If you're looking for a higher level language to do this, then consider MATLAB (or its open source work-alike, Octave.) Actually, Maple can do this kind of computation too, but it can be quite hard to force Maple into using floating point arithmetic rather than its normal symbolic computation mode.

Answer by Brian Borchers for Nonstandard Hessian approximations in Gauss-Newton

2012-09-06T05:04:04Z

You could certainly try this, but you'd have to do a lot of careful analysis to derive any convergence results. Among other things to consider:

Your objective E(x) is more likely to have local minima due to the nonconvexity introduced by the logarithms. This could also happen with the sum of squares objective, but in practice it's uncommon for reasonably well behaved f(x).
Your approximate Hessian will typically be fully dense, and depending on the size of $x$, this might make the solution of the equations impractical.

You haven't said anything about how you're computing $J$ or how you might be able to compute the second derivatives of $f$.

You haven't said anything about how large your vector of parameters $x$ is. If it is large, then rather than attempting this, I'd suggest using a limited memory BFGS method to avoid storing the dense $H$.

If it's small, then depending on the difficulty of computing the second derivatives of $f$, I'd probably use the full Hessian and implement Newton's method rather than doing all of the work to show that this Gauss-Newton like method had good theoretical properties. If the second derivatives are hard to compute, then I'd just use a conventional BFGS quasi-Newton method.

Answer by Brian Borchers for Levenberg-Marquadt near the minima for non-zero-residual problems

2012-09-04T00:33:19Z

The technique of multiplying the diagonal by $(1+\lambda)$ fails when you've got a 0 on the diagonal of $J^{T}J$. Having a zero on the diagonal of $J^{T}J$ can happen when you're far away from the optimal solution (in which case it tells you nothing about what might be true at optimality), or it can happen at an optimal solution (where the singularity of $J^{T}J$ can be a suggestion that the optimal solution is not unique.)

As Thomas stated in his answer, a common approach is to use the QR factorization to solve the least squares problem at each iteration and zero out small diagonal elements of $R$.

If you're stuck using the Cholesky factorization (e.g. if your problem is large and sparse and the Cholesky factorization can reasonably be computed but the QR factorization can't be computed), then you can avoid this problem with zeros on the diagonal by adding a small multiple of the identity to $J^{T}J$ rather than by multiplying the diagonal elements by $(1+\lambda)$.

You should also look carefully at the scaling of your problem- bad scaling frequently leads to poor convergence in practice.

Answer by Brian Borchers for Cascading minimization problems

2012-08-11T17:35:29Z

You can solve the first LP, obtain the optimal value, and then add a constraint that the first objective has that optimal value. Then you can add your second set of constraints to the original set of constraints and solve that LP. e.g. if your original LP is

$\max c^{T}x $

subject to

$Ax=b$

$x \geq 0$

Suppose the optimal objective value is $p^{*}$. Now suppose you want to solve a second problem

$\max d^{T}x $

subject to

$ Fx=g$

$ x \in U$

This can be formulated as

$\max d^{T}x $

subject to

$Ax=b$

$c^{T}x=p^{*}$

$Fx=g$

$x \ge 0$

Most simplex based LP codes can efficiently reoptimize after adding the additional constraints.

Answer by Brian Borchers for linear programming with OR restrictions

2012-08-11T14:26:01Z

Your additional constraints make the feasible region for your problem nonconvex, and thus it cannot be represented as a linear programming problem.

If you can obtain an upper bound $M_{i}$ on the maximum value of $x_{i}$, then there is a standard approach to these problems in which the problem is formulated as a 0-1 mixed integer linear programming problem.

First, we introduce 0-1 variables $y_{i}$, and add the constraints

$x_{i} \geq k_{i}y_{i}$

If $y_{i}$ takes on a 1 value, then $x_{i}$ is forced to be greater than or equal to $k_{i}$. If $y_{i}$ is 0, then this constraint is vacuous.

Next, we add the constraints

$x_{i} \leq M_{i}y_{i}$

If $y_{i}$ is 0, this forces $x_{i}=0$. If $y_{i}=1$, then this constraint does nothing.

This combination of constraints means that either $x_{i}=0$ and $y_{i}=0$ or $x_{i}\geq k_{i}$ and $y_{i}=1$.

There are many approaches to solving the resulting 0-1 mixed integer linear programming including branch and bound methods and cutting plane algorithms. In practice, the most powerful methods (implemented in closed source commercial codes such as IBM's CPLEX as well as a number of open source noncommercial software packages) combine these two general approaches into a "branch and cut" approach.

It is also possible to directly implement these "semi-continuous variables" within a branch and cut algorithm, and this approach does not require an upper bound $M_{i}$. This feature is available for example in IBM's CPLEX package.

Answer by Brian Borchers for Characterising semi-definite positiveness on vectors with non-negative entries

2012-07-24T17:10:28Z

Your cone $C$ is the cone of copositive matrices. The dual of C is the cone of compeltely positive matrices. See e.g.

http://mathworld.wolfram.com/CopositiveMatrix.html

Answer by Brian Borchers for Ease of calculation of norm

2012-07-16T23:45:09Z

I'm assuming that by $\| A^{1/2}b-z \|$, you're referring to the 2-norm and that by $A^{1/2}$, you're referring to the unique symmetric matrix square root of $A$.

If you can precompute $A^{1/2}z$, then you can quickly compute $\|A^{1/2}b-z\|$.

$\| A^{1/2}b-z \|_{2}^{2}=(A^{1/2}b-z)^{T}(A^{1/2}b-z)$

$\| A^{1/2}b-z \|_{2}^{2}=b^{T}Ab-2b^{T}A^{1/2}z + z^{T}z$

$\| A^{1/2}b-z \|_{2}=\sqrt{b^{T}Ab-2b^{T}A^{1/2}z + z^{T}z}$

Note that in many cases, the Cholesky factorization of $A$ can be used in place of the symmetric matrix square root.

Answer by Brian Borchers for SDP Algorithms/ maximally complementary solutions

2012-07-09T16:03:35Z

The standard approach is to embed the original SDP in a self dual formulation that has strictly feasible primal/dual solutions, solve the self dual formulation and then reach conclusions about the original problem from the solution of the self-dual problem.
See for example:

http://citeseerx.ist.psu.edu/viewdoc/download?doi=10.1.1.37.679&rep=rep1&type=pdf Infeasible-start semidefinite programming algorithms via self-dual embeddings

Unfortunately, I believe that this will only tell you when the problem has no strictly complementary primal-dual solutions and not necessarily give you a maximally complementary solution.

Answer by Brian Borchers for On Random Vectors and Eigenvectors of Symmetric Matrices

2012-07-03T01:58:57Z

In this case, it's easy enough to compute the probability exactly in terms of the incomplete beta function.

Let $v$ be the fixed unit vector, and $r$ be the random unit vector, uniformly distributed over the $n$-dimensional hypersphere.

$P(| r^{T} v | < \alpha) = 1-P(| r^{T}v | \geq \alpha) $

$P(| r^{T} v | < \alpha) = 1-2P(r^{T}v \geq \alpha) $

Let $\theta$ be the angle between $r$ and $v$. Then $\cos \theta=r^{T}v$. Now, $r^{T}v \geq \alpha$ only if $r$ is in the spherical cap of the unit hypersphere centered around $v$, containing vectors within an angle $\cos^{-1}\alpha$ from $v$. The height of this spherical cap is $h=1-\alpha$.

The probability that $r^{T}v \geq \alpha$ is then given by the ratio of the surface area of this circular cap to the surface area of the hypersphere. Using standard formulas for the surface area of the hypersphere and of the spherical cap, we get that

$P(| r^{T}v | < \alpha)=1-I_{2h-h^{2}}(\frac{n-1}{2},\frac{1}{2})$

A quick test with a Monte Carlo simulation in MATLAB verifies this formula.

Answer by Brian Borchers for Efficient Algorithm For Projection Onto A Convex Set

2012-06-26T17:25:57Z

You haven't told us anything about the size of your problem instances. How many terms are there in the sum of norms? What is $n$?

Your problem is an example of a "sum of norms" optimization problem. Searching with Google Scholar will lead you to published research on this class of problems.

CVX is using a standard approach for solving this problem by reformulating it as a second order cone programming (SOCP) problem and then using a primal-dual interior point method (SeDuMi's or SDPT3) to solve the resulting SOCP. For small problem instances, this should be a very robust and reasonably fast approach to solving the problem, but there are faster primal-dual codes for SOCP available (both CPLEX and MOSEK can be used to solve SOCP's.) You could ask CVX to extract the SOCP problem and export it in SeDuMi format and then try to use another solver on the SOCP.

You might also look at first order methods to solve the SOCP- For example, I believe that TFOCS could be used to solve the problem.

Answer by Brian Borchers for Moore-Penrose bound question

2012-06-09T13:26:50Z

If $Ax=b$ has a unique solution $x^{*}$, then $x_{m}=x^{*}$.

If $Ax=b$ has infinitely many solutions, then $x_{m}$ will be one of these solutions. In particular, it will be the solution with the smallest two-norm. There will be other solutions $x$ such $\| x-x_{m} \|$ is arbitrarily large.

If $Ax=b$ has no solutions then it isn't clear at all what you mean by referring to a specific $x$.

Are you asking about what happens in the presence of noise in $A$ or the right hand side vector $b$? In these cases the condition number of $A$ can be used to bound the effect of the noise on the pseudoinverse solution. Unfortunately, if $\mbox{cond}(A)=\infty$, there is no bound.

Answer by Brian Borchers for Convolutive noise removal

2012-06-03T04:16:31Z

A lot depends on $\hat{\eta}(\xi)$. When you convolve this with $\hat{u}$ and $l(\xi)$, you will lose the sparsity if $\hat{\eta}(\xi)$ has broad support.

Just how much do you know about the spectrum of $\eta(t)$? If you know it well enough, you can deconvolve it before trying your basis pursuit approach. If you don't know it very well, then you're in deep trouble.

Answer by Brian Borchers for My overdetermined linear system gives both bad and good estimates. Why ?

2012-05-30T19:16:49Z

In using least squares, you normally want the residuals from the various equations to be indpendent of each other. If your $q_{i}$ vectors are imprecise measurements, then this will introduce correlation between the residuals in the pair of equations invovling $q_{i}$.

When you say that you get a "good fit" using the homoegneous system, how well does the solution satisfy the original nonhomogeneous system of equations?

Answer by Brian Borchers for Multiple Linear Regression Estimation without full recalc

2012-05-28T16:39:55Z

There's a very large literature on updating solutions to least squares problems as new data are added. The naive formula $\hat{\beta}=(X^{T}X)^{-1}X^{T}y$ can be problematic in practice because of numerical issues with ill-conditioning. The QR factorization is typically a much better choice in practice. There are lots of papers on updating the QR factorization as data are added. See for example:

http://eprints.ma.man.ac.uk/1192/01/covered/MIMS_ep2008_111.pdf

Answer by Brian Borchers for A question about normalization in the Fourier transform

2012-05-28T03:51:15Z

There are many variations of the definition of the Fourier transform and its inverse that are all essentially equivalent. For example, the factor of $2\pi$ often appears as $1/(2\pi)$ in front of the integral in one of the transform pair, or as $1/\sqrt{2\pi}$ in front of both integrals in the transform pair.

Researchers in different fields have often adopted a particular convention. If you want to publish in a journal where such a convention is well established, then you will probably find that it is easier to get your paper accepted for publication if you follow the convention used by people working in that field.

Answer by Brian Borchers for Nonlinearly constrained optimization (quadratic)

2012-05-21T18:47:37Z

The real issue here is the constraint

$\sum_{i} x_{i}1_{x_{i}>a} < b $

whose left hand side has horrible discontinuities.

Rather than using a solver designed for problems with continuous variables, you should formulate this as 0-1 mixed integer nonlinear programming problem, with binary decision variables $z_{i}$ that are 0 when $x_{i} \leq a$ and 1 when $x_{i}>a$. You can then use branch and bound to solve the problem.

Answer by Brian Borchers for Continuity of Lexicographic Minimum Solution of a parametrized LP problem

2012-05-19T21:28:13Z

You've chosen a somewhat limited form of "parametric LP" here. It's more common to allow $F$ to vary as $F=F_{0}+tF_{1}$. For that more general variety of parametric LP, the conjecture is false- it's easy to construct a counterexample. For example, consider

$\min tx_{1}+x_{2}$

$x_{1}+x_{2}=1$

$x \geq 0$

What happens as t varies from $t<1$ to $t=1$ to $t>1$?

You could also consider parametric changes in $A$- the conjecture is also false for that case.

Returning to the problem as stated:

It's easy to construct examples where the LP becomes infeasible at some values of $t$. You'll need to exclude such values of $t$.

Since this posting looks like it could well be a homework exercise, I'll simply outline an approach that you might use to prove this. First note that the basic feasible solutions (BFS's) for your parametric system of inequalities are (except at values of $t$ where a BFS vanishes or first appears) continuous functions of $t$. Also note that for each value of $t$, the set of optimal solutions is the convex hull of these BFS's. Finally, show that $x_{L}(t)$ is always at a BFS. Now, go back and think about what happens at those values of $t$ where BFS's vanish or appear.

Answer by Brian Borchers for Solving for an operator by minimization

2012-05-04T19:06:25Z

You haven't said so, but I'm assuming that $\psi$ and $\phi$ are vectors. These could more generally be functions in some function space, and you would typically discretize those functions to work with vectors in your numerical computations.

I'm also assuming that you have one or more $\psi$, $\phi$ examples containing measured values or values that have been computed by some numerical model. Call these examples $\psi_{1}$, $\phi_{1}$, $\psi_{2}$, $\phi_{2}$, $\ldots$, $\psi_{n}$,$\phi_{n}$.

There may well be no exact solution to $O(x,y)\psi=\phi$ that works for all of these examples, so a typical approach would be to turn this into a nonlinear least squares problem of the form

$\min_{x,y} \sum_{k=1}^{n} \| O(x,y)\psi_{k}-\phi_{k} \|_{2}^{2} $

Your data and parameters are complex valued, but these are most easily dealt with by splitting the real and imaginary parts to produce a nonlinear least squares problems involving only real parameters and measured values.

Once you've reduced the problem to a conventional real nonlinear least squares problem, you can use standard methods to solve it. The Levenberg-Marquadt method is most commonly used in practice. In the LM method, you can use finite difference approximations to get the required derivatives- most software for doing this has the ability to do this finite differencing for you.

Answer by Brian Borchers for Explicit formula for Cholesky factorization in a special case

2012-04-25T01:23:59Z

The matrix $\alpha J$ is a rank one matrix, so there are simple update/downdate formulas for computing the Choleksy factorization of $Q+sI-\alpha J$ if you start with the factorization of $Q+sI$.

I'm not aware of any update formulas that get you from the Cholesky factorization of $Q$ to a Cholesky factorization of $Q+sI$.

Answer by Brian Borchers for Solving a particular nonlinear system of equalities

2012-04-24T00:18:17Z

Your problem might be small enough that it is within the range of polynomial optimization techniques based on SDP relaxations of sums of squares problem. This has been implemented in software packages such as GLOPTIPOLY. See

http://homepages.laas.fr/henrion/software/gloptipoly3/

Answer by Brian Borchers for Ways to convert a Positive Semi-Definite (PSD) matrix -> Positive Definite matrix

2012-04-19T14:00:59Z

The choice that you make can result in a huge difference in the solution. Neither method is particularly good and both can be quite unstable.

For example, suppose that

$M=[1\;0\;0\; 0;\; 0 \; 1 \; 0 \; 0; \; 0 \; 0 \; 0 \; 0;\; 0 \; 0 \; 0 \; 0]$

and

$N=[1;\; 1; \; 1; \; 1]$.

What do you want $L$ to be in this situation? Why?

Now, suppose that you try using

$\delta = 1 \times 10^{-14}$.

Then your approximation to $M^{1/2}$ is

$M^{1/2}=[1\;0\;0\; 0;\; 0 \; 1 \; 0 \; 0; \; 0 \; 0 \; 1 \times 10^{-7} \; 0;\; 0 \; 0 \; 0 \; 1 \times 10^{-7}]$

and your approximation to $M^{-1/2}$ is

$M^{-1/2}=[1\;0\;0\; 0;\; 0 \; 1 \; 0 \; 0; \; 0 \; 0 \; 1 \times 10^{7} \; 0;\; 0 \; 0 \; 0 \; 1 \times 10^{7}]$.

In this case, the least squares solution is

$L=[1; \; 1; \; 1 \times 10^{14} \; 1 \times 10^{14}]$.

In the same example, if you use the pseudoinverse, then you'd get

$L=[1; \; 1; \; 0 \; 0]$.

Using the pseudoinverse effectively eliminated from consideration any projection of $N$ onto the null space of $M$. That might be what you want to do.

However, if $M$ is ill-conditioned, then it becomes practically impossible to distinguish small eigenvalues of $M$ from $0$ eigenvalues. As a result the computation using the pseudoinverse also becomes extremely unstable and dependent on the tolerance used in computing the pseudoinverse. This isn't good either.

Continuing the example, suppose that

$M=[1\;0\;0\; 0;\; 0 \; 1 \; 0 \; 0; \; 0 \; 0 \; 1 \times 10^{-14} \; 0;\; 0 \; 0 \; 0 \; 1 \times 10^{-14}]$.

The pseudoinverse solution is then

$L=[1; \; 1; \; 1 \times 10^{14};\; 1 \times 10^{14}]$

Now, suppose that due to a rounding tolerance in computing the pseudoinverse, you compute

$\mbox{pinv}(M)=[1\;0\;0\; 0;\; 0 \; 1 \; 0 \; 0; \; 0 \; 0 \; 0 \; 0;\; 0 \; 0 \; 0 \; 0]$

Then you'd get

$L=[1; \; 1; \; 0 \; 0]$

again.

You really need to back up a bit and explain to us why you're trying to solve this particular least squares problem in a situation where $M$ might be singular. It's likely that you need to deal with this possible ill-conditioning at a higher level by using some sort of regularization approach that is based on the particular aspects of your underlying problem.

Answer by Brian Borchers for Sherali-Adams relaxation

2012-04-18T16:19:20Z

It isn't clear from your posting whether you're trying to understand:

Why the inequalities generated by the Sherali-Adams procedure are valid?

Why the procedure is complete in the sense that after enough iterations you arrive at the convex hull of the integer solutions of the original integer linear programming problem?

I suppose that you might also be interested in the question of how this can be used in practice.

You should be aware that Sherali-Adams is just one of several "lift and project" schemes based on LP or SDP relaxations that provide a ladder of relaxations of an integer programming problem from the simple LP relaxation up to a relaxation that has only integer solutions as its extreme points.

You might find the following paper by Monique Laurent to be useful in understanding this stuff:

http://citeseerx.ist.psu.edu/viewdoc/summary?doi=10.1.1.18.2521

Answer by Brian Borchers for sparsity of QR decomposition

2012-04-16T13:33:12Z

Unfortunately you can't. With any orthogonal factorization (e.g. QR, LQ, or SVD) you have the problem that because some of the columns of the orthogonal matrix have to span a particular subspace, and because the remaining columns have to form an orthogonal basis for the complement to this subspace, and because these spaces can be completely arbitrary, the orthogonal matrix won't be sparse unless you happen to be very lucky.

There are "sparse QR" methods that effectively represent Q as a product of Givens rotations rather than storing Q explicitly. The Givens rotations are transformations that are extremely sparse/structured matrices.

Answer by Brian Borchers for Solving a non-convex quadratically constrained quadratic program

2012-04-12T15:49:34Z

Branch-and-Bound can be used to solve problems like this, but it's an exponential time algorithm that is not practical for large instances of the problem. What is $n$, the dimension of the $x$ vector? How big is your $m$?

Comment by Brian Borchers

2013-05-23T00:25:15Z

The general result is that if you have a quadratically constrained quadratic programming problem with only one non-convex constraint, then you can formulate it exactly as an SDP. That's why I believe that this relaxation will be exact in this case. However, I haven't looked into this for several years, so the details escape me at the moment. If I get a chance I'll try to sit down and work through the details.

Comment by Brian Borchers

2013-05-21T14:20:39Z

See appendix B of "Convex Optimization" by Boyd and Vandenberghe.

Comment by Brian Borchers

2013-05-21T04:02:31Z

On second thought, although approach I sketched out is sufficient to show that the SDP is a relaxation of the original problem, it's not enough to show that the relaxation is exact (I forgot to account for the $L^{T}CL$ term in the constraint. I'm afraid that you need the machinery of the S-lemma to show that the relaxation is exact.

Comment by Brian Borchers

2013-05-21T03:35:35Z

As I recall, you can also get to this formulation using the S-lemma.

Comment by Brian Borchers

2013-05-21T03:23:07Z

There are several software packages for solving such an SDP, including SDPA (written in C++) and CSDP (in C)

Comment by Brian Borchers

2013-03-05T04:31:53Z

Your notation seems a bit confused (and perhaps suggests why you're unable to establish this result.) You haven't explained what $a_{0}$, $a_{1}$, $...$ are.

Comment by Brian Borchers

2013-02-03T04:19:21Z

The poster has also put this same question up on scicomp.stackexchange.com, where she's somewhat more likely to get a useful answer, provided that she can clarify the question.

Comment by Brian Borchers

2012-12-26T03:44:33Z

You can think of the $x_{i}^{2}$ as nonnegative weights that sum to one. This opens up the whole world of the generalized mean inequality.

Comment by Brian Borchers

2012-12-16T14:25:24Z

If the matrix isn't sparse, and the cost of getting individual matrix entries is large compared to the cost of accessing an element of a matrix stored in conventional dense matrix form, then iterative methods are going to be horribly slow in practice.

Comment by Brian Borchers

2012-12-16T14:22:25Z

Let me clarify what I meant here- "being able to get an arbitrary element M(i,j) at little cost" isn't very useful. If you don't know where the nonzero elements are in the matrix, then you have to check every single one to find the nonzeros. If you do happen to know where the nonzero elements are, and you can compute them quickly, then you could use this as a way to do matrix vector multiplications in an iterative method.

Comment by Brian Borchers

2012-10-11T03:42:57Z

See ams.org/samplings/feature-column/… for a discussion of the Stern-Brocot tree which has been used to design optimal gear trains. Tables of optimal gear trains based on this method have been around for a long time.

Comment by Brian Borchers

2012-09-07T13:08:16Z

On the other hand, if this really an unconstrained problem, I'd simply use BFGS. The computations inside the BFGS method won't be time consuming for a small problem like this, so the big concern will be the cost of function and derivative computations. You'll probably want to experiment with finite difference derivatives vs. derivatives by automatic differentiation vs. analytical formulas for the derivatives to see what works best (both in terms of accuracy and speed.) You can also play with the convergence criteria- perhaps a relatively imprecise solution will be adequate for your work.

Comment by Brian Borchers

2012-09-07T04:18:26Z

Now I'm even more confused about your problem. You say that x lives "in a non-euclidean manifold", so it appears that you don't have an unconstrained optimization problem. In that case, you'd have a constrained nonlinear optimization problem and you'd need to consider algorithms for constrained problems rather than methods for unconstrained optimization. A good (but somewhat old) textbook that deals with mathematics and many practical aspects of optimization is "Practical Optimization" by Gill, Murray, and Wright. You'd find lots of good advice in that book on how to approach this.

Comment by Brian Borchers

2012-09-05T18:01:45Z

@Alex; The $J^{T}J$ approximation to the hessian of $c(x)$ works well when $f_{i}(x)-y_{i}$ is small, but if the residuals are large then the approximation can degrade. This is because the second order terms that go into the Hessian are dropped when we approximate it with $J^{T}J$, and those terms have $f_{i}(x)-y_{i}$ factors.

Comment by Brian Borchers

2012-09-05T17:59:00Z

This is confusing- are you minimizing a sum of squares or are you using some other objective (such as the Huber measure)? If you're not minimizing a sum of squares than LM isn't appropriate in the first place. Assuming that you are doing nonlinear least squares, and assuming that you have a reasonable assessment of the uncertainty in the measurements (e.g. $y_{i}$ is known with uncertainty $\sigma_{i}$, then you should normalize by $c(x)=\sum ((f_{i}(x)-y_{i})/\sigma_{i})^{2}$. You should also scale the parameters $x$ so that they're all of similar magnitude.