11

Given code that computes a function $f(x)$, automatic differentiation tools produce a code that can compute $f(x)$ and its derivatives at the same time. Solving a differential equation is an entirely different problem and AD doesn't solve differential equations (although AD tools are sometimes useful in connection with PDE constrained optimization.) AD ...


7

Julia has a whole ecosystem for generating sparsity patterns and doing sparse automatic differentiation in a way that mixes with scientific computing and machine learning (or scientific machine learning). Tools like SparseDiffTools.jl, ModelingToolkit.jl, and SparsityDetection.jl will do things like: Automatically find sparsity patterns from code Generate ...


6

Optim.jl from Julia will work with the number types that you give it, so if you make it use BigFloats then it'll do that. Local derivative based, derivative-free, global, and integrates with automatic differentiation. From Julia, it's just: using Optim rosenbrock(x) = (1.0 - x[1])^2 + 100.0 * (x[2] - x[1]^2)^2 result = optimize(rosenbrock, big.(zeros(2)), ...


6

First off, let me make a subtle distinction in the implementation of level set methods in topology optimization. In the literature, you will see this method implemented using shape derivatives as in here or using material derivatives (through the use of a Heaviside function) as in here. In my experience, using shape derivatives works better, but material ...


5

The disciplined convex rules used by CVX don't know how to deal with ratios or products of variables or functions of variables, so this can't be written directly in CVX. However, it can be written using $\mbox{rel_entr}$ in CVX. To see this, $y \; \mbox{entr}(x/y)= -y (x/y) \log(x/y) = -x \log(x/y) = -\mbox{rel_entr}(x,y)$ where $\mbox{rel_entr}(x,y)=x \...


5

Because you say the losses are convex, I will presume that all $c_i \ge 0$, which means that max is used in a convex fashion. Given that, this problem can be formulated as a Linear Programming problem (LP). Define additional optimization variables, $y_i$. Replace $f_i(x_i)$ with $c_iy_i$, and add the constraints $y_i \ge x_i, y_i \ge 0$. The result is an LP. ...


5

This is just an NLP (Non-Linear Programming) model. You can rewrite it as: $$\begin{align}\min \>& Z \\ & \sum_i w_i = 1 \\ & \sum_i f_i(w_i\cdot Z)\ge k\cdot Z \\ & w_i \in [0,1] \end{align}$$ Getting rid of a division is always a good idea. If we can assume $Z\gt0$, then a slightly different formulation can look like: $$\begin{align}\...


4

The general idea that you have of learning an easy to compute model from results of your detailed simulation model and then optimizing the easy to compute model is long-established. The easy to compute model is typically called a surrogate model or a response surface model. Once the surrogate is available, you can use conventional optimization techniques ...


4

Let $f(\omega)$ be your power spectrum. Then maybe something like $$ \frac{\|f\|_{L^\infty}\|f\|_{L^0}}{\|f\|_{L^1}} = \frac{\mathrm{max}_{\omega\in\Omega} f(\omega)\cdot|\omega_{max}-\omega_{min}|}{\int_\Omega|f(\omega)|d\omega}. $$ I know that $L^0$ isn't really good notation but I think it is useful for presenting this. This quantity is minimized by ...


4

This topic has been discussed at some length on Cross Validated (aka stats.stackexchange) and Reddit: Why is Newton's method not widely used in machine learning? (see in particular Nick Alger's answer) Why use gradient descent with neural networks? L-BFGS and neural nets Why second order SGD convergence methods are unpopular for deep learning? How does the ...


4

If your objective function is noisy, then it makes sense to use stochastic algorithms. I would take a look at James Spall's SPSA algorithm and variations. There are also algorithms that update (an approximation of) the Hessian. All of these algorithms take into account that whatever gradient you compute may be noisy, and consequently don't just blindly ...


4

If you have a good LP solver, then the linear programming approach often works well. You don’t want to implement your own simplex or interior point code for LP though. A specialized variant of the simplex method due to Barrodale and Roberts is a popular approach to this problem, but it takes some effort to implement this efficiently and you might be better ...


3

Despite my comment, I think you can find $\tilde{D}$ that contains the noise term as well. You have this equation: $$-M^{T} \tilde{D} M \phi(t) = -M^{T} D M \phi(t) + W(t)$$ Where $W(t)$ is the noise term vector. So: $$-M^{T} (\tilde{D}-D) M \phi(t) = W(t)$$ Take $\mathcal{D} = \tilde{D} - D$. Let's expand this equation: $$(M \phi(t))_{i} = \sum_{j=1}^{...


3

The determinant is the product of all eigenvalues $\lambda_i(A)$, so its logarithm is the sum of the logarithms of the eigenvalues. As a consequence, you can write the objective function as follows: $$ \min_A \quad-\sum_i \log\lambda_i(A). $$ Now, assume you had an eigenvalue that tried to go to zero, then its logarithm would go to negative infinity and so ...


3

Following up on my comment on the original question, I have finally managed to construct a counter example that shows that the statement is not in fact correct. Define the positive part of a function, $$ [x]^+ = \begin{cases}x & \text{if $x\ge 0$} \\ 0 & \text{otherwise.}\end{cases} $$ The let $$ \sigma_1(y) = 1+ \left[\tfrac 14 - |y-1|\right]^+ $...


2

3D local feature descriptors for shapes are very well studied. Typically, people tend to represent the input as a set of points (point clouds) and try to characterize the local neighborhoods with lower dimensional signatures a.k.a. descriptors. The traditional descriptors, analogous to their 2D counterparts, involve some kind of histogram-ing such as SHOT, ...


2

You should use a modeling language so that your code is independent of the underlying solver. cvxpy is a good choice. When I rewrite your model in cvxpy: #!/usr/bin/env python3 import cvxpy as cp import numpy as np from numpy.random import normal as randn Sample = 10 H = randn(size=(4,2,Sample))+1j*randn(size=(4,2,Sample)) h = randn(size=(4,1))+1j*randn(...


2

I guess a problem with this approach is that it also "sparsifies" the gradient. If you look at the objective function: \begin{align} \Phi(x) = \lvert|(Ax)\odot(Ax)-b |\rvert^2+\lambda\lvert|x|\rvert_1 \end{align} If one uses a proximal gradient method, we have \begin{align} x_{k+1} = \mathcal{P}\left(x_k-\alpha A^T\left(\left((Ax_k)\odot(Ax_k)-b\...


2

Is it what you need? (λ is unknown) # INPUT N,K = 5,3 h=[1.1,1.2,1.7,0.5,0.3] Xk = cp.Variable()# unknown Lambda Xn = [cp.Variable() for i in range(N-K)] constraints = [] for x in Xn: constraints.append(cp.abs(x)<=Xk) X = [Xk]*K+Xn obj = cp.Minimize(cp.sum([(h[i]-X[i])**2 for i in range(N)])) prob = cp.Problem(obj, ...


2

What about this: Step 1: For each edge $P_iP_j$, look at the two triangles that share it: $P_iP_jP_k$ and $P_lP_jP_i$. Step 2: Compute the counter-clockwise normal on each of the two triangle: $$\vec{n}_{ijk} = \frac{\vec{P_iP_j} \times \vec{P_jP_k}}{\|\vec{P_iP_j} \times \vec{P_jP_k}\|} \,\,\, \text{ and }\,\,\, \vec{n}_{lji} = \frac{\vec{P_lP_j} \times ...


2

My usual answer is "don't use fitness proportionate selection". If you want to use it though, you kind of have to enter the world of tuning things to get the level of selection pressure you want. You could, instead of $1/r(x)$, do $(k+1)/(k+r(x))$ for some problem-specific value of $k$. That'll scale things to some degree. You could apply some non-linear ...


2

You will want to look at the chapter on "Hessian modification" methods in the excellent book "Numerical Optimization" by Nocedal and Wright. You will find that the Levenberg-Marquardt method is probably what you are looking for, given that you can't easily compute eigenvalues of large sparse matrices.


2

So I went ahead and just modified the scipy optimize.py source to do what I wanted. I assume this is not best practices, but... moving on. So I used the inspect module to figure out where the corresponding functions are. So for the Powell optimization routine, I found: def _minimize_powell on line 2509 of scipy/optimize/optimize.py. Then I went to line ...


2

50 is a lot of parameters. You could try doing a basic first order sensitivity analysis to determine whether you can drop any of these. Using Bayesian Optimization to minimize a cost function is one way of dealing with the problem you've encountered. But remember that your standard L2 norm might have counterintuitive behaviours in high dimensions (see On the ...


2

Think of an object that moves along a spiral -- say, an electron moving in a uniform magnetic field. Its positions are not periodic (it never comes back to its original place) but its velocities are. Looking at periodicity in the velocities therefore tells us something we don't know by just looking at the positions.


2

If a certain velocity component is periodic with period $\tau$ that means that the corresponding coordinate, as a function of time, is a sum of a linear function and a periodic function with the same period. To prove it, let $\dot{x}(t)$ be periodic, $\dot{x}(t+\tau)=\dot{x}(t)$, and the integral over the period $\int_{t}^{t+\tau} \dot{x} dt = I$, where $I$ ...


2

Thanks to Wolfgang Bangerth for pointing out that this can be rewritten as a linear problem. My reformulation would be: $$\min_{x\in\mathbb{R}^N}\|Ax-b\|_1 \implies$$ $$\min_{(x,y)\in\mathbb{R}^{N+M}}\sum_{i=1}^N y_i, $$ $$Ax - y \leq b$$ $$-(Ax+y) \leq -b$$ $$y\geq0$$


2

You should have a look at Satisfiability Modulo Theories or SMT for short. A huge number of problems can be thought of as instances of SMT for a particular theory. For example, correctly designing certain types of integrated circuits can be phrased as an SMT problem. The problem you're describing fits under the theory of quantifier-free linear integer ...


1

CVXPY's norm atom won't accept a raw Python list as an argument; you need to pass it a CVXPY expression. Stack the list of scalars into a vector using the hstack atom, like so: constraints = [cp.norm( cp.hstack([ y_hat[col] - cp.trace( np.transpose((B_hat_star[:,col][:,np.newaxis]*np.sqrt(L)*C_hat[col,:])) @ X) for col in ...


1

The condition for finding a root that corresponds to positive definiteness of the Hessian in optimization, is that the function grows strictly monotonically. But this is not a useful condition. That is because even in optimization, positive definiteness of the Hessian does not actually guarantee convergence of the unmodified Newton method. What is important ...


Only top voted, non community-wiki answers of a minimum length are eligible