Recent Questions - Computational Science Stack Exchange most recent 30 from scicomp.stackexchange.com 2022-07-06T18:43:00Z https://scicomp.stackexchange.com/feeds https://creativecommons.org/licenses/by-sa/4.0/rdf https://scicomp.stackexchange.com/q/41577 0 Which preconditioners make Richardson iteration convergent? Yaroslav Bulatov https://scicomp.stackexchange.com/users/18786 2022-07-06T17:38:55Z 2022-07-06T18:02:11Z <p>Suppose we solve an <span class="math-container">$m\times n$</span> full-rank system of equations <span class="math-container">$Ax=b$</span> by iterating the following for a small enough <span class="math-container">$\mu&gt;0$</span></p> <p><span class="math-container">$$x=x+\mu B(b-Ax)$$</span></p> <p>Is there a nice description of kinds of <span class="math-container">$B$</span> which make this iteration convergent? For instance, for <span class="math-container">$n=1$</span>, <span class="math-container">$A,B$</span> can be viewed as vectors, and this iteration converges iff their dot product is positive. What about <span class="math-container">$n&gt;1$</span>?</p> <ol> <li><span class="math-container">$B=A^\dagger$</span> is (trust-region) Newton's method</li> <li><span class="math-container">$B=A^T$</span> is gradient descent on least-squares objective</li> <li><span class="math-container">$B=I$</span> is modified Richardson iteration, only works for positive-definite <span class="math-container">$A$</span></li> </ol> <p>Curious if there's a family of methods between 2. and 3. -- converges for any full-rank <span class="math-container">$A$</span>, but without relying too much on <span class="math-container">$A'y$</span> (much more expensive than <span class="math-container">$Ax$</span> in my application)</p> https://scicomp.stackexchange.com/q/41576 0 Centered finite volume scheme for an advective term on an unstructured/irregular/non-uniform grid nicholaswogan https://scicomp.stackexchange.com/users/38057 2022-07-06T15:38:47Z 2022-07-06T16:53:56Z <p>Consider the continuity equation</p> <p><span class="math-container">$$\frac{\partial u}{\partial t} + \frac{\partial \Phi}{\partial x} = 0$$</span> <span class="math-container">$$\Phi = au + b\frac{\partial u}{\partial x}$$</span></p> <p>Suppose I want to solve the above using centered finite volumes. Below, <span class="math-container">$U_j \approx \overline{u}_j$</span> where <span class="math-container">$\overline{u}_j$</span> is the spatial average value of <span class="math-container">$u$</span> in a finite volume. <span class="math-container">$\Phi_{i+1/2}$</span> and <span class="math-container">$\Phi_{i-1/2}$</span> are the fluxes at the edges of the finite volume, and <span class="math-container">$\Delta x_i$</span> is the thickness of the the volume.</p> <p><span class="math-container">$$\frac{\partial U_j}{\partial t} = -\frac{\Phi_{i+1/2} - \Phi_{i-1/2}}{\Delta x_i}$$</span></p> <p>For a structured grid, I can use the following centered approximation for the fluxes</p> <p><span class="math-container">$$\Phi_{i+1/2} = a_{j+1/2}\frac{U_{j+1} + U_{j}}{2} + b_{j+1/2}\frac{U_{j+1} - U_{j}}{\frac{1}{2}\Delta x_{j+1} + \frac{1}{2}\Delta x_{j}}$$</span></p> <p>For the advective term above, I have approximated <span class="math-container">$U_{j+1/2}$</span> as</p> <p><span class="math-container">$$U_{j+1/2} = \frac{U_{j+1} + U_{j}}{2}$$</span></p> <p><strong>My question is what is a more general centered approximation of <span class="math-container">$U_{j+1/2}$</span> for an unstructured grid?</strong> One possibility is to linearly interpolate to <span class="math-container">$U_{j+1/2}$</span> between the center of the finite volumes containing <span class="math-container">$U_{j+1}$</span> and <span class="math-container">$U_{j}$</span>. The result is</p> <p><span class="math-container">$$U_{j+1/2} = \frac{U_{j+1} \Delta x_j + U_{j} \Delta x_{j+1}}{\Delta x_{j+1} + \Delta x_{j}} \tag{1}$$</span></p> <p>A reference equation in a textbook would be appreciated.</p> <p>Thanks!</p> https://scicomp.stackexchange.com/q/41574 0 Using absolute error as the cost function Natasha https://scicomp.stackexchange.com/users/29087 2022-07-06T14:54:45Z 2022-07-06T14:54:45Z <p>This is related to my previous post <a href="https://scicomp.stackexchange.com/questions/41552/minimize-distance-between-curves">Minimize distance between curves</a>.</p> <p>I have a dataset with values of multiple curves. An example plot is shown below. I want to scale the curves (move up/down) so that all curves overlap.</p> <p>The following is a sample dataset which includes that data points corresponding to 5 curves and coordinate inputs below</p> <pre><code>scale = 1.5; x1 = [0,4,6,10,15,20]*scale; y1 = [18,17.5,13,12,8,10]; x2 = [0,10.5,28]*scale; y2= [18.2,10.6,10.3]; x3 = [0,4,6,10,15,20]*scale; y3 = [18,13,15,12,11,9.6]; x4 = [9,17,28]*scale; y4 = [5,5.5,7]; x5 = [1,10,20]*scale; y5 = [3,0.8,2]; plot(x1,y1, '*-', x2, y2, '*-', x3, y3, '*-', x4, y4, '*-', x5, y5, '*-') </code></pre> <p><a href="https://i.stack.imgur.com/WTTld.png" rel="nofollow noreferrer"><img src="https://i.stack.imgur.com/WTTld.png" alt="enter image description here" /></a></p> <p>To scale the curves, I need to find the scale factor for each curve.</p> <p>In the answer posted in the previous post (<a href="https://scicomp.stackexchange.com/a/41555/29087">https://scicomp.stackexchange.com/a/41555/29087</a>) the following has been suggested.</p> <p><span class="math-container">$$\min_{a} \left(\sum_i ( a \cdot f_i -g_i)^2\right)$$</span> which leads to a single equation <span class="math-container">$$0 = \partial_a \left(\sum_i ( a \cdot f_i -g_i)^2\right) = 2 \sum_i f_i (a \cdot f_i -g_i)$$</span> and thus <span class="math-container">$$a = \frac{\sum f_i g_i}{\sum f_i^2}$$</span></p> <p>In the above answer, the cost function is the sum of absolute error, and differentiating this gives <code>a</code>.</p> <p>Excuse me for the naive question,</p> <p>I would like to understand how the cost function has to be defined if we want to minimize absolute error instead of squared error and how to estimate the scale factor.</p> <p>Suggestions will be really appreciated.</p> https://scicomp.stackexchange.com/q/41570 1 Finding weighted average of curves Natasha https://scicomp.stackexchange.com/users/29087 2022-07-05T17:22:48Z 2022-07-06T07:34:28Z <p>This is related to my previous post <a href="https://scicomp.stackexchange.com/questions/41552/minimize-distance-between-curves">here</a></p> <p>I have a dataset with values of multiple curves. An example plot is shown below. I want to scale the curves (move up/down) so that all curves overlap.</p> <p>The following is a sample dataset which includes that data points corresponding to 5 curves and coordinate inputs below</p> <pre><code>scale = 1.5; x1 = [0,4,6,10,15,20]*scale; y1 = [18,17.5,13,12,8,10]; x2 = [0,10.5,28]*scale; y2= [18.2,10.6,10.3]; x3 = [0,4,6,10,15,20]*scale; y3 = [18,13,15,12,11,9.6]; x4 = [9,17,28]*scale; y4 = [5,5.5,7]; x5 = [1,10,20]*scale; y5 = [3,0.8,2]; plot(x1,y1, '*-', x2, y2, '*-', x3, y3, '*-', x4, y4, '*-', x5, y5, '*-') </code></pre> <p><a href="https://i.stack.imgur.com/WTTld.png" rel="nofollow noreferrer"><img src="https://i.stack.imgur.com/WTTld.png" alt="enter image description here" /></a></p> <p>To scale the curves, I need to find the scale factor by defining a target curve. I'm not sure of the ways in which the target curve can be defined. Would it be a good approach to compute the weighted average? Since the x scale is different for each curve, I am not sure how to define an average/ target curve.</p> <p>Suggestions will be really appreciated.</p> https://scicomp.stackexchange.com/q/41568 0 What are the Exact Rules for Significant Figures, Precision, and Uncertainty? CoastCity Lapse 00crashtest https://scicomp.stackexchange.com/users/43486 2022-07-05T06:37:21Z 2022-07-06T15:07:17Z <p>In the physical sciences (which are physics, chemistry, astronomy, materials science, etc.), we learned that the uncertainty is +/- the smallest unit (which is 1) of the last significant figure if the uncertainty is not given in a recording of data. So, if we have a digital measuring device that measures to the nearest millimeter, has a manufacturer's stated uncertainty of +/- 1 mm, and gives a reading of 914 mm, then it will obviously be recorded as just &quot;914 mm&quot;.</p> <p>However, does the true value actually lie somewhere between exactly 913 mm and exactly 915 mm, or may it stray outside even those numbers if higher precision is used? For example, if go down to the micrometer, is the uncertainty actually +/- 999 μm or +/- 1,499 μm according to the rules of significant figures? If we measure the same sample using a micrometer, is the reading guaranteed to be somewhere between 913,001 microns and 914,999 microns, or is it instead only guaranteed to be somewhere between 912,501 microns and 915,499 microns, respectively?</p> https://scicomp.stackexchange.com/q/41565 3 My Complex Matrix SVD is Correct according to rule A = USV' but Wrong according to Matlab or any linear algebra library Emir evcil https://scicomp.stackexchange.com/users/43495 2022-07-04T20:21:59Z 2022-07-05T18:43:26Z <p>I am working on Singular Value Decomposition for complex matrices. I implemented One Sided Jacobi algorithm. It gives exactly the same result as the svd function in Matlab for the real matrices. However, it does not give the same result in complex matrices. It confirms the formula A = UEV* however I am not sure whether this result is correct or not. Any idea why the result is different? Can I say that the result I obtained is absolutely true if it meets the relation A = USV*?</p> <p>Do you have any suggestions about algorithms that work stably in complex matrices? Do you know which SVD algorithms using on libraries such as Matlab, Scipy, LAPACK are?</p> <p>My and Numpy outputs:</p> <p><strong>A:</strong></p> <p>[(0.563612363551724+0.4358699131559025j), (0.2806793692533188+0.9930311421704819j), (0.5202215785270062+0.4481106844815471j)] [(0.0011901478455322856+0.6075877153377014j), (0.9997530325963236+0.464449485996173j), (0.8705886268046512+0.6556403407040208j)]<br /> [(0.13985961755470633+0.2579120041586833j), (0.09334365980890702+0.0107473125009665j), (0.9198536550293809+0.6672661113638818j)]</p> <p><strong>U According to my result:</strong></p> <p>[(0.32263888225141446+0.49316146595502114j), (0.24956005697294448+0.42711101441813026j), (0.5117539922520571-0.39255845198236494j)]<br /> [(0.5230826883717128+0.4215011626318952j), (0.28920165671330744-0.01699244263936318j), (-0.7012897455688379+0.17167125279675066j)]<br /> [(0.3775020380130938+0.24271259233586367j), (-0.5162076548934359-0.6363169611154663j), (0.24545327028669411+0.04980039176599202j)]</p> <p><strong>S according to my result :</strong></p> <p>2.228961805416005 0.7603266403486695 0.6863428620403935</p> <p><strong>V according to my result:</strong></p> <p>[(0.4048700775539+0.15820544389995128j), (0.4186462523569478-0.1933392387422724j)(0.850657804557617+0.019852985890792853j)] [(0.5434747415164244+0.03930807376934979j), (0.6282038478106724-0.0013258064898205727j), (-0.559364893463487+0.03670246154675491j)]<br /> [(0.7614928175839254-0.11216861617097869j), (-0.6938072848356582-0.11786215892893791j), (-0.019427233541509478-0.1874673042557197j)]</p> <p><strong>U according to numpy and Matlab:</strong></p> <p>[[-0.16336677-0.53547895j 0.66551039-0.11000351j -0.47255786-0.09087257j] [-0.37024417-0.60367276j -0.3400191 -0.14303676j 0.44455477-0.40595198j] [-0.22386647-0.36736338j -0.24426998+0.59097786j -0.10864474+0.62785478j]]</p> <p><strong>S according to numpy and MATLAB:</strong></p> <p>[2.26594193 0.91577243 0.47884249]</p> <p><strong>V according to numpy and MATLAB:</strong></p> <p>[[-0.36133177+0.j -0.55295939+0.19926312j -0.65938003+0.29864268j] [ 0.39102224+0.j -0.37701504+0.67597439j 0.08382929-0.49090992j] [-0.84648743+0.j 0.06187986+0.2271988j 0.32018676-0.35424718j]]</p> <p><strong>USV' :(conjugate transpose of V) according to my results -&gt;</strong></p> <p>[(0.5636123635517245+0.43586991315590246j), (0.2806793692533189+0.9930311421704825j), (0.520221578527007+0.44811068448154784j)]<br /> [(0.0011901478455318415+0.6075877153377023j), (0.9997530325963246+0.4644494859961733j), (0.8705886268046524+0.6556403407040218j)]<br /> [(0.13985961755470683+0.25791200415868387j), (0.09334365980890735+0.010747312500966673j), (0.9198536550293819+0.6672661113638825j)]</p> <p>A = USV' is confirmed but i could not understand why results are different from ready libraries.</p> <p>Since my code in python is a bit complicated, I leave the implementation of the same method in matlab here.</p> <pre><code>function [U,S,V]=jacobi_svd1(A) % Floating Point function equivalent to MATLAB function svd(A) implemented % using 1-sided Jacobi algorithm [m,n]=size(A); % Get size of matrix A U = A; % Assign U as A V=eye(n); % Assign V as identity matrix of size n count=5; % Number of sweeps %% while(count&gt;=1) for i = 1:n-1 for j = i+1:n a = norm(U(:,i),2); % Calculate norm of ith column b = norm(U(:,j),2); % Calculate norm of jit column % Assure the singular values will appear in decreasing order in S % swap columns i and j of U and V if a &lt; b temp(:,j) = U(:,j); U(:,j) = U(:,i); U(:,i) = temp(:,j); temp1(:,j) = V(:,j); V(:,j) = V(:,i); V(:,i) = temp1(:,j); end %% % Compute submatrix of U'U x=0; y=0; w=0; for k=1:m x=x+(U(k,i))^2; end for k=1:m y=y+(U(k,j))^2; end for k=1:m w=w+(U(k,i))*(U(k,j)); end %% % Compute the Jacobi rotation that diagonalizes the % submatrix if w ~= 0 alpha=(y-x)/(2*w); if alpha&gt;=0 t = 1/(abs(alpha)+sqrt(1+alpha^2)); else t = -(1/(abs(alpha)+sqrt(1+alpha^2))); end c=1/sqrt(1+t^2); s=c*t; else c=1; s=0; end %% % update columns i and j of U T = U(:,i); U(:,i)=c*T-s*U(:,j); U(:,j)=s*T+c*U(:,j); % update matrix V of right singular vectors T = V(:,i); V(:,i)=c*T-s*V(:,j); V(:,j)=s*T+c*V(:,j); end end count = count - 1; end %% %singular values are the norms of the columns of U for j=1:n singvals(j)=norm(U(:,j),2); if (singvals(j) &gt; 0) U(:,j) = U(:,j)/singvals(j); end end S=diag(singvals); % Arrange singular values along the diagonal of S end </code></pre> <p>Code reference : <a href="https://github.com/jaysmitjadhav/Singular-Value-Decomposition/blob/master/SVD/jacobi_svd1.m" rel="nofollow noreferrer">https://github.com/jaysmitjadhav/Singular-Value-Decomposition/blob/master/SVD/jacobi_svd1.m</a></p> https://scicomp.stackexchange.com/q/41564 2 Is there a permutation used in sparse QR factorizations that better locates small elements on the diagonal of R? wyer33 https://scicomp.stackexchange.com/users/10094 2022-07-04T18:54:54Z 2022-07-04T18:54:54Z <p>Is there a permutation used in sparse QR factorizations that better locates small elements on the diagonal of R to the end of the diagonal? As an example, consider the following snippet of MATLAB code:</p> <pre><code>A = [1 1 1; ... -1 -1 1; ... 0 0 -1; ... 0 0 1; ... -1 0 0; ... 0 -1 0]; s = [1e-8;1e-8;1e-8;1;1e-8;1]; B = sparse([A diag(s)]); [Q,R,p] = qr(sparse(B)'); pp = [5,6,1,2,4,3]; [QQ,RR] = qr(sparse(B(pp,:)')); </code></pre> <p>This matrix is similar to what may occur in an optimization problem with small slack variables. If we look at the absolute value of the diagonal elements of <code>R</code>, we see a small element in position 5:</p> <p><a href="https://i.stack.imgur.com/bv9KO.png" rel="nofollow noreferrer"><img src="https://i.stack.imgur.com/bv9KO.png" alt="Diagonal elements of R" /></a></p> <p>I would prefer this element to be at the end of the diagonal of <code>R</code>. Using the permutation <code>pp</code>, this can be achieved and the the diagonal of <code>RR</code> is seen in this image:</p> <p><a href="https://i.stack.imgur.com/NwRdR.png" rel="nofollow noreferrer"><img src="https://i.stack.imgur.com/NwRdR.png" alt="enter image description here" /></a></p> <p>In this case, the sparsity between <code>R</code> and <code>RR</code> is the same. Ideally, I'd like a diagonal sorted by absolute value similar to what is achieved in a dense QR factorization. In that case, it's easy to get a rough idea of the effective rank of the original matrix and truncate the matrix <code>R</code> effectively for the applications that I care about. At the same time, I understand that this would dramatically impact the sparsity of <code>R</code> and for large matrices this is not tractable. In truth, I don't really need a fully sorted diagonal, but some kind of permutation that mostly puts the small elements at the end. At the moment, I have some very large matrices and may get some very small elements in the first few dozen diagonal elements, which is undesirable.</p> <p>As some additional information, if it helps, I understand that MATLAB uses <code>spqr</code> as its underlying QR factorization algorithm and that the parameter <code>spqrtol</code> can be set in <code>spparms</code> to adjust the drop tolerance for small elements. Alternatively, spqr could just be used directly. In this case, I'd contend the issue is how to effectively set this tolerance if the issue is small diagonal elements compared to other diagonal elements.</p> <p>Thanks for any insights.</p> https://scicomp.stackexchange.com/q/41563 0 Upper bound for conditional relative entropy Botao Hao https://scicomp.stackexchange.com/users/32559 2022-07-04T18:50:52Z 2022-07-04T18:50:52Z <p>I am wondering if the following inequality could hold up to some constants: <span class="math-container">$$\mathbb E_X[D_{KL}(P_{Y|X}, Q_{Y|X})]\leq c D_{KL}(P_{X, Y}, Q_{X,Y}),$$</span> where <span class="math-container">$c&gt;0$</span> is some constant, <span class="math-container">$P$</span> and <span class="math-container">$Q$</span> are two probability measures.</p> <p>The difficulty now is that KL does not have triangle inequality. If anyone knows whether this is true, I appreciate it a lot!</p> https://scicomp.stackexchange.com/q/41562 1 I need MATLAB or Mathematica code to solve this integral over limit 0- λ tolulope ojuola https://scicomp.stackexchange.com/users/0 2022-07-04T08:28:59Z 2022-07-04T15:17:19Z <p>I want the code to integrate equation(1) or (2) over the limits using mathematical or mathlab to get equation (3) as the answer of <span class="math-container">$Z$</span> vibrational partition function, giving the following additional information <span class="math-container">\begin{aligned} Z_{v i b}(\beta) &amp;=\sum_{n=0}^{\lambda} e^{-\beta E_{n}}, \quad \beta=\frac{1}{k T}, \\ \lambda &amp;=-c+\sqrt{A} \pm \sqrt{A-B} \end{aligned}</span>, <span class="math-container">$$Z_{v i b}(\beta)=\sum_{n=0}^{\lambda} e^{\frac{\beta a^{2} \hbar^{2}}{4 m}(B-2 A)+\frac{\beta a^{2} \hbar^{2} B^{2}}{8 m(n+c)^{2}}+\frac{\beta a^{2} \hbar^{2}}{8 m}(n+c)^{2}}$$</span> In the classical limit, the summation turns to integral <span class="math-container">$$\begin{array}{l} Z_{v i b}(\beta)=\int_{0}^{\lambda} e^{\frac{\beta a^{2} \hbar^{2}}{4 m}(B-2 A)+\frac{\beta a^{2} \hbar^{2} B^{2}}{8 m(n+c)^{2}}+\frac{\beta a^{2} \hbar^{2}}{8 m}(n+c)^{2}} d n,............(1) \\ =\int_{c}^{\lambda+c} e^{\frac{\beta a^{2} \hbar^{2}}{4 m }(B-2 A)+\frac{\beta a^{2} \hbar^{2} B^{2}}{8 m \rho^{2}}+\frac{\beta a^{2} \hbar^{2}}{8 m \prime} \rho^{2}} d \rho, \rho=c+n ...........(2)\\ =e^{\beta \Lambda_{3}} \sqrt{\pi}\left[\frac{\left[1-\operatorname{Erf}\left(\Lambda_{1} \sqrt{\beta}\right)\right]-e^{\Lambda_{4} \beta}\left[1-\operatorname{Erf}\left(\Lambda_{2} \sqrt{\beta}\right)\right]}{\Lambda_{5} \sqrt{\beta}}\right]...........(3) \end{array}$$</span></p> <p>where we have also introduced the following parameters for mathematical simplicity, <span class="math-container">$$\begin{array}{l} \Lambda_{1}=\frac{1}{\rho} \sqrt{-\frac{a^{2} \hbar^{2} B^{2}}{8 m}}-\rho \sqrt{-\frac{a^{2} \hbar^{2}}{8 m}}, \quad \Lambda_{2}=\frac{1}{\rho} \sqrt{-\frac{a^{2} \hbar^{2} B^{2}}{8 m}}+\rho \sqrt{-\frac{a^{2} \hbar^{2}}{8 m}} \\ \Lambda_{3}=-\frac{a^{2} \hbar^{2} A}{2 m}, \Lambda_{4}=\frac{a^{2} \hbar^{2} B}{2 m}, \Lambda_{5}=\sqrt{-2 a^{2} \hbar^{2}}, \quad c \leq \rho \leq c+\lambda \end{array}$$</span></p> https://scicomp.stackexchange.com/q/41561 -1 how to solve multisource minimum total completion time problem in network Dustin https://scicomp.stackexchange.com/users/43492 2022-07-04T07:03:30Z 2022-07-04T07:03:30Z <p>Assume that there is a graph G_t(V,E_t) at every time slot t, V is the nodes in the graph and E_t is the edges that link two nodes, every edge e has a capacity C_e which means how much data can pass through the edge at the same time. Some nodes S={s_1,s_2,...,s_n} are the sources in V, who have data to transmit to one destination dst. The data volume that needs to be transferred to the dst of source s_i is d_i. Then we want to find out the paths for every source at every time slot and the corresponding data rate on each path to minimize the time that all sources finish transferring all the data. Note that every source can use multiple paths at the same time but no more than \phi.</p> <p>What is the best way to solve this problem? what about the complexity of this problem, is it an np-hard problem?</p> https://scicomp.stackexchange.com/q/41560 0 Computing the overlapping area of curves Natasha https://scicomp.stackexchange.com/users/29087 2022-07-04T04:58:47Z 2022-07-04T05:11:19Z <p>I have a set of datapoints that correspond to 5 curves. Sample dataset and MATLAB inputs:</p> <pre><code>scale = 1.5; x1 = [0,4,6,10,15,20]*scale; y1 = [18,17.5,13,12,8,10]; x2 = [0,10.5,28]*scale; y2= [18.2,10.6,10.3]; x3 = [0,4,6,10,15,20]*scale; y3 = [18,13,15,12,11,9.6]; x4 = [9,17,28]*scale; y4 = [5,5.5,7]; x5 = [1,10,20]*scale; y5 = [3,0.8,2]; plot(x1,y1, '*-', x2, y2, '*-', x3, y3, '*-', x4, y4, '*-', x5, y5, '*-') </code></pre> <p>I want to find the overlapping area under these curves.</p> <p><a href="https://i.stack.imgur.com/mEBEg.png" rel="nofollow noreferrer"><img src="https://i.stack.imgur.com/mEBEg.png" alt="enter image description here" /></a> Suggestions on the functions that could be used to find the intersecting area and the x-limits of the intersecting area will be really helpful.</p> https://scicomp.stackexchange.com/q/41558 4 How does one proceed to solve (big) underdetermined or overdetermined systems of linear equations "nowadays"? Sen90 https://scicomp.stackexchange.com/users/43488 2022-07-03T17:33:29Z 2022-07-05T20:02:33Z <p>In my numerical linear algebra class we mentioned this problem briefly and according to some other lectures on the internet especially in data driven environments one mostly has to deal with such over- or underdetermined systems. Nevertheless in our lecture we mostly deal with direct methods (mostly matrix decompositions) and iterative methods (CG, GMRES, MINRES, ...) for different types of square! matrices.</p> <p>Are there similar methods for under- or overdetermined systems?</p> https://scicomp.stackexchange.com/q/41554 -1 How to initialize Eigen C++ parameters within for-loop? Shakil https://scicomp.stackexchange.com/users/43482 2022-07-02T10:39:26Z 2022-07-02T20:17:53Z <p>I am new to Eigen C++ programming. I am trying to create an Rcpp function to call from R. The function takes a list of matrices (Xlst), and two lists of vectors (ylst and smwlst). In each for-loop iteration, the function should perform standardization of the columns of a matrix and then multiply with the corresponding vectors. Please see the code below.</p> <p>It appears that I am not initializing the parameters correctly, or I should initialize the parameters outside of the for-loop, which I don't know how to do because the number of rows are not the same for all matrices. Could anyone help with this?</p> <p>Rcpp function:</p> <pre><code>List solveR_obj_rLog2_cov(List ylst, List Xlst, List smwlst, int P, int H) { int P, H, h, j, it=0; Eigen::MatrixXd BetaH = Eigen::MatrixXd::Zero(H, P); for(h = 0; h &lt;= H; ++h) { Eigen::VectorXd yh = ylst[h]; Eigen::VectorXd lamh = smwlst[h]; Eigen::MatrixXd Xh = Xlst[h]; Nh = yh.size(); double numh, dnmh; Eigen::VectorXd yhc(Nh), Xhm(P), Xhsd(P); Eigen::MatrixXd Xhs(Nh,P); yhm = yh.mean(); yhc = yh.array() - yhm; for(j = 0; j &lt; P; ++j) { Xhm(j) = Xh.col(j).mean(); Xh.col(j) = Xh.col(j).array() - Xhm(j); Xhsd(j) = sqrt(Xh.col(j).squaredNorm()/Nh); Xhs.col(j) = Xh.col(j)/Xhsd(j); } for(j = 0; j &lt; P; ++j) { numh = (lamh*yhc.dot(Xhs.col(j))).sum()/Nh; dnmh = (lamh*Xhs.col(j).dot(Xhs.col(j))).sum()/Nh; BetaH.coeffRef(h,j) = numh/dnmh; } } return BetaH; } </code></pre> <p>Data from the R:</p> <pre><code>yl &lt;- list(allData20$y, allData21$y, allData22$y) smwl &lt;- list(allData20$w, allData21$w, allData22$w) Xl &lt;- list(as.matrix(allData20[,2:5]), as.matrix(allData21[,2:5]), as.matrix(allData22[,2:5])) </code></pre> <p>Application in the R:</p> <pre><code> c_res &lt;- solveR_obj_rLog2_cov(ylst = yl, Xlst = Xl, smwlst = smwl, P = 4, H = 3) </code></pre> https://scicomp.stackexchange.com/q/41552 4 Minimize distance between curves Natasha https://scicomp.stackexchange.com/users/29087 2022-07-01T17:03:06Z 2022-07-06T11:43:35Z <p>I have a dataset with values of multiple curves. An example plot is shown below. I want to shift the curves (up/down) so that all curves overlap. This would mean the data points in each curve is scaled up/down by a factor. I am not sure how to frame this as a mathematical problem (to minimize the vertical distance between each pair of curves) and determine the scaling factor for each curve. I tried to start by computing the pair-wise distance matrix but I am not sure what to do next.</p> <p>Suggestions will be really appreciated.</p> <p><a href="https://i.stack.imgur.com/Vref3.png" rel="nofollow noreferrer"><img src="https://i.stack.imgur.com/Vref3.png" alt="enter image description here" /></a></p> <p>EDIT:</p> <p>The following is a sample dataset which includes that datapoints corresponding to 5 curves and coordinate inputs are below</p> <pre><code>scale = 1.5; x1 = [0,4,6,10,15,20]*scale; y1 = [18,17.5,13,12,8,10]; x2 = [0,10.5,28]*scale; y2= [18.2,10.6,10.3]; x3 = [0,4,6,10,15,20]*scale; y3 = [18,13,15,12,11,9.6]; x4 = [9,17,28]*scale; y4 = [5,5.5,7]; x5 = [1,10,20]*scale; y5 = [3,0.8,2]; plot(x1,y1, '*-', x2, y2, '*-', x3, y3, '*-', x4, y4, '*-', x5, y5, '*-') </code></pre> <p>And the plot looks like below. I tried to do a piecewise interpolation using a linear approximation to find the function values. Since the curves don't have a common x, I am not sure how the target function has to be selected.</p> <p><a href="https://i.stack.imgur.com/2yvOYm.png" rel="nofollow noreferrer"><img src="https://i.stack.imgur.com/2yvOYm.png" alt="enter image description here" /></a></p> <p>Suggestions will be really appreciated. Also, I only want to scale the datapoints <span class="math-container">$$\tilde f_i(a) = a \cdot f_i$$</span> and do not want to shift it along the x-dimension (I'm not sure if <code>b</code> <span class="math-container">$\tilde f_i(a,b) = a \cdot f_i$</span> + b shifts the function value along the x-dimension; may be I am wrong and <code>b</code> shifts the function value in the y-direction). And I also prefer to define a tolerance window/interval <span class="math-container">$\epsilon$</span> (user-defined) above and below the target function which will allow scaling the other curves such that they don't completely overlap on the target curve.</p> <p>EDIT2: I would like to ask for a suggestion. If I want to scale <span class="math-container">$f_i$</span> (move it a bit up) and also scale <span class="math-container">$g_i$</span> (move it a bit down) would it make sense if the objective function is expressed as <span class="math-container">$(a⋅fi − gi)^2 + (b.gi - fi)^2$</span>. Basically, I want to try and move both the curves towards each other instead of moving one curve towards the other.</p> <p>EDIT3: @Vladislav Gladkikh</p> <p>I like this idea, this preserves the shape of all <span class="math-container">$x_k$</span>'s i.e the trend of the <span class="math-container">$x_k$</span>'s before and after shifting doesn't change and I find this great. Could you please add a pictorial explanation of the third point?</p> <p>I think overall we do, <span class="math-container">$f(x_k) - \bar{f(x_k)} + \bar{g(x_k)}$</span> Here, f is the non-main curve and g is the main curve (i.e y values of the curves at the corresponding x values), and <span class="math-container">$\bar{f(x_k)}$</span> is the average.</p> <blockquote> <p>Then shift all curves toward the main curve so that their averages coincide within the corresponding domain segments.</p> </blockquote> <p>I think I am a bit confused here since we are subtracting/adding the mean the <span class="math-container">$\bar{f(x_k)}$</span>'s/ <span class="math-container">$\bar{g(x_k)}$</span>'s computed at the <span class="math-container">$x_k$</span>'s in the entire domain. I am not sure if my interpretation of segment is wrong, please correct me if I am wrong.</p> <p>Could you please clarify if we refer to a segment as the domain between any 2 points that lie on the curve xk, yk here. Could you please illustrate this in a figure, if possible?</p> https://scicomp.stackexchange.com/q/41549 3 Finite difference on matrices PC1 https://scicomp.stackexchange.com/users/41372 2022-07-01T04:03:48Z 2022-07-01T04:03:48Z <p>This question relates closely to <a href="https://scicomp.stackexchange.com/q/41541/41372">another question</a> I asked earlier this week.</p> <p>Let's say that I want to find a set of matrices <span class="math-container">$S_i$</span>, with <span class="math-container">$0\leq i\leq N$</span>, that minimizes some objective function with numerous terms. I can easily compute the gradient of the objective function, with regard to each matrix. This provides me a system of equations where each matrix is linked to its immediate neighbours.</p> <p>For example, in my case, we can write <span class="math-container">$S_i=F_i(S_{i-1}, S_{i+1})$</span>, where <span class="math-container">$F_i$</span> is some function. We also have some conditions at the boundaries so <span class="math-container">$S_0$</span> and <span class="math-container">$S_N$</span> are fully determined.</p> <p>In some sense, this is very similar to a 1D finite difference problem, except that the <span class="math-container">$S_i$</span> are matrices, not humble scalars.</p> <p>For scalars, the problem is normally simple enough to write, for example as a linear problem with one big (and eventually sparse) matrix. There are many ways to solve this type of problem. However, for matrices, it can get quite challenging as matrices do not commute in general so the function <span class="math-container">$F$</span> can be more complex to linearize and write as a large matrix, especially if it also involves the inverse of some matrices <span class="math-container">$S_i$</span>.</p> <p>Is there any resource on how to solve these kind of problems? Or at least to reformulate it in a way that is more consistent with how people look at this type of problems? For now the best I can do is to solve each matrix from its neighbours over and over, so it takes a long time to eventually converge.</p> https://scicomp.stackexchange.com/q/41548 0 Solving PDE with a non-linear constraint in MATLAB penghao zhang https://scicomp.stackexchange.com/users/43476 2022-07-01T01:39:03Z 2022-07-01T12:59:47Z <p>I am trying to solve a DAE with a non-linear constraint. The governing equations have the following form.</p> <p><a href="https://i.stack.imgur.com/3LuJh.png" rel="nofollow noreferrer"><img src="https://i.stack.imgur.com/3LuJh.png" alt="enter image description here" /></a></p> <p>The second equation is a constraint and it must be satisfied everywhere. Is there a way to use the ode solver in MATLAB to solve this type of system? If so, for example, ode15s, how could we couple the constraint into the mass matrix?</p> https://scicomp.stackexchange.com/q/41547 3 How can I check mass conservation when solving the advection equation using an upwind scheme? nicholaswogan https://scicomp.stackexchange.com/users/38057 2022-06-30T19:08:04Z 2022-07-05T16:31:55Z <p>My question is how to keep track of the &quot;mass&quot; being advected out of a model domain, for the 1-D advection equation, and an upwind differencing scheme. Following is the background</p> <p>Consider the continuity equation</p> <p><span class="math-container">$$\frac{\partial u}{\partial t} + \frac{\partial \Phi}{\partial x} = 0$$</span></p> <p>Suppose the flux <span class="math-container">$\Phi$</span> is given by an advective term</p> <p><span class="math-container">$$\Phi = a u$$</span></p> <p>Suppose <span class="math-container">$a = 0.01x$</span>, and we are considering only <span class="math-container">$x&gt;0$</span>. Therefore, <span class="math-container">$a &gt; 0$</span>, so the wave is traveling in the positive <span class="math-container">$x$</span> direction. I can discretize the spatial derivative with a first order upwind scheme. <span class="math-container">$j$</span> is grid cell centers, and <span class="math-container">$j-1/2$</span> is the left edge of the grid cell.</p> <p><span class="math-container">$$\frac{\partial \Phi}{\partial x}|_j \approx \frac{\Phi_{j} - \Phi_{j-1/2}}{\frac{1}{2}\Delta x}$$</span></p> <p><span class="math-container">$$= \frac{a_{j} u_j - a_{j-1/2} u_{j-1/2}}{\frac{1}{2}\Delta x}$$</span></p> <p>I replace <span class="math-container">$u_{j-1/2}$</span> with the averge of the grid cell neighbors (<span class="math-container">$\frac{u_{j-1} + u_{j}}{2}$</span>)</p> <p><span class="math-container">$$= \frac{a_j u_j - a_{j-1/2} \frac{u_{j-1} + u_{j}}{2}}{\frac{1}{2}\Delta x}$$</span></p> <p>So, using the method of lines, I've now turned the PDE into a system of ODEs,</p> <p><span class="math-container">$$\frac{\partial u_j}{\partial t} = - \frac{2 a_j u_j - a_{j-1/2} u_j - a_{j-1/2} u_{j-1}}{\Delta x} \tag{1}$$</span></p> <p>Now we can consider the left boundary (<span class="math-container">$j = 1$</span>)</p> <p><span class="math-container">$$\frac{\partial \Phi}{\partial x}|_1 = \frac{\Phi_{1} - \Phi_{1/2}}{\frac{1}{2}\Delta x}$$</span></p> <p>Here, <span class="math-container">$\Phi_{1/2} = \Phi_\text{in}$</span> is the &quot;flux&quot; entering the left side of the model domain. The left boundary is then given by</p> <p><span class="math-container">$$\frac{\partial u_1}{\partial t} = \frac{-2 a_1 u_1}{\Delta x} + \frac{2 \Phi_\text{in}}{\Delta x} \tag{2}$$</span></p> <p>Equations (1) and (2) describe a system of ODEs which can be evolved forward in time. I have done this for an initial condition of <span class="math-container">$u_j = 0$</span>, and <span class="math-container">$\Phi = 0.1$</span> for the domain <span class="math-container">$x = [0,100]$</span>. The model reaches the following steady state:</p> <p><a href="https://i.stack.imgur.com/8Mt1p.png" rel="nofollow noreferrer"><img src="https://i.stack.imgur.com/8Mt1p.png" alt="results" /></a></p> <p>Mass is being advected out of the right side of the model domain. How do I figure out this mass flux? The mass flux out should be the same as <span class="math-container">$\Phi_\text{in} = 0.1$</span>, for flux conservation. Any other relevant tips for understanding/checking mass conservation in these scenarios is appreciated.</p> <p>I would have thought <span class="math-container">$\Phi_\text{out} = u_na_n$</span>, where <span class="math-container">$n$</span> is the right most grid cell, but this is <span class="math-container">$\approx 0.105$</span>. So, mass is not quite being conserved, or I'm not computing <span class="math-container">$\Phi_\text{out}$</span> properly.</p> <p>Thanks!</p> https://scicomp.stackexchange.com/q/41541 5 Recurrence relation for matrices PC1 https://scicomp.stackexchange.com/users/41372 2022-06-29T18:48:04Z 2022-07-01T09:22:20Z <p>I have matrices (<span class="math-container">$S_0$</span> thought <span class="math-container">$S_N$</span>) and I have a recurrence relation that link successive matrices together. <span class="math-container">$$S_i + S_i(aS_{i-1}^{-1})S_i=C_i+aS_{i+1}$$</span></p> <p>We can assume for this problem that <span class="math-container">$S_0=S_N=I$</span>, the identity matrix. The <span class="math-container">$C_i$</span> in this equation are all known beforehand. All matrices <span class="math-container">$S_i$</span> are definite positive and <span class="math-container">$a&gt;0$</span>. The matrices <span class="math-container">$C_i$</span> are positive semi-definite, with some eigenvalues equal to 0 but at least one eigenvalue is positive.</p> <p>If I look at the problem with <span class="math-container">$1\times1$</span> matrices, the analytic solution for <span class="math-container">$S_t$</span> in terms of <span class="math-container">$S_{t-1}$</span> and <span class="math-container">$S_{t+1}$</span> is quite straightforward to write and there is always exactly one positive root. I would like to see if I can find an iterative scheme to compute the matrices <span class="math-container">$S_i$</span>, even if it requires many iterations and it takes forever...</p> <p><strong>Edit:</strong> The simplified question is to compute <span class="math-container">$X$</span> from <span class="math-container">$X+XAX=B$</span>, where <span class="math-container">$A$</span> and <span class="math-container">$B$</span> are known. All matrices here are symmetrical and definite-positive.</p> <p><strong>Edit 2:</strong> Looking at the equation <span class="math-container">$X+XAX=B$</span>, it is very similar to a <a href="https://en.wikipedia.org/wiki/Quadratic_equation" rel="nofollow noreferrer">quadratic equation</a>. So I thought that maybe we can solve it in a similar way.</p> <p>We can factor the equation as: <span class="math-container">\begin{align}XAX+X&amp;=B\\ \left(A^\frac12X+\frac12A^{-\frac12}\right)^T\left(A^\frac12X+\frac12A^{-\frac12}\right)&amp;=\frac14A^{-1}+B\\ &amp;=\frac{A^{-1}}4\left(I+4AB\right)\end{align}</span> The matrix <span class="math-container">$A$</span> is positive-definite so <span class="math-container">$A^\frac12$</span>, <span class="math-container">$A^{-\frac12}$</span>, and <span class="math-container">$A^{-1}$</span> all exist and also are positive-definite. That means that <span class="math-container">$I+4AB$</span> is also positive-definite. So we can &quot;extract the square root&quot; and rewrite: <span class="math-container">\begin{align}\left(A^\frac12X+\frac12A^{-\frac12}\right)&amp;=\frac{A^{-\frac12}}2\left(I+4AB\right)^\frac12\\ X+\frac12A^{-1}&amp;=\frac{A^{-1}}2\left(I+4AB\right)^\frac12\\ X&amp;=\frac{A^{-1}}2\left[-I+\left(I+4AB\right)^\frac12\right] \end{align}</span></p> <p>This is basically the same approach as the one used to solve a quadratic equation, except that I solve a matrix equation. Also, it's clear that the positive root only is used, so we can retrieve a matrix <span class="math-container">$X$</span> that is also definite-positive.</p> <p><strong>Updated question:</strong> Is this approach valid? I know I shouldn't do an updated question but I will copy it in an answer only if it's valid...</p> https://scicomp.stackexchange.com/q/41540 5 Software for Feasibility Problems Dan Doe https://scicomp.stackexchange.com/users/36438 2022-06-29T10:20:31Z 2022-06-30T08:48:44Z <p>I face a feasibility problem of type</p> <p><span class="math-container">$$c_i(\boldsymbol x) \leq 0, i = 1, \dots, \mathcal{I} \\ c_e(\boldsymbol x) = 0, e = 1, \dots, \mathcal{E}$$</span></p> <p>where <span class="math-container">$\mathcal{I} + \mathcal{E} \gg \text{dim}(\boldsymbol x) \sim \mathcal{O} (10^1)$</span>.</p> <p>Currently, I solve this with <a href="https://coin-or.github.io/Ipopt/" rel="nofollow noreferrer">Ipopt</a> but since it takes quite some time (many iterations reqiuired) I thought about looking for some special feasible point solvers.</p> <p>The &quot;largest&quot; (3) collection of feasible point solvers I could find online are mentioned in section 5.3 of <a href="http://www2.math.uni-wuppertal.de/opt/preprints/prep2010/amna_opap_10_03.pdf" rel="nofollow noreferrer">this paper</a>. Unfortunately, there seems to be no implementation of the therein developed algorithm (<code>EFNES</code>) online. The <a href="https://www.galahad.rl.ac.uk/doc/filtrane.pdf" rel="nofollow noreferrer"><code>FILTRANE</code></a> framework of the <a href="https://github.com/ralna/GALAHAD" rel="nofollow noreferrer">Galahad package</a> might be quite effective but has a nasty Fortran API which looks even for the example problem quite difficult.</p> <p>EDIT: The link in the original <code>TRESNEI</code> paper is dead, but I found <a href="https://flore.unifi.it/handle/2158/1271224" rel="nofollow noreferrer">this working one</a>.</p> <p>Any suggestions besides the usual suspects from nonlinear optimization, ideally free to use within academia?</p> https://scicomp.stackexchange.com/q/41533 0 Solving differential equations with fast oscillations using odeint danial https://scicomp.stackexchange.com/users/43454 2022-06-27T16:45:58Z 2022-06-29T15:02:02Z <p>I have wrote this code to solve an equation , I know the behavior of this function has very rapid oscillations, when I RUN it gives bogus values for some &quot;m[x]&quot; and some &quot;t&quot;'s, with this error:</p> <blockquote> <p>C:\Users\dani\anaconda3\lib\site-packages\scipy\integrate\odepack.py:247: ODEintWarning: Excess work done on this call (perhaps wrong Dfun type). Run with full_output = 1 to get quantitative information. warnings.warn(warning_msg, ODEintWarning)</p> </blockquote> <pre><code>import scipy as sio import numpy as np import mpmath as mp import scipy.integrate as spi import matplotlib.pyplot as plt import time initial_value=np.logspace(24,27,100) t=np.logspace(0,6,100) m=np.logspace(0,6,100) start_time=time.perf_counter() phi_m={} phi_m_prime={} phi=[] phi_prime=[] j=0 i=np.pi*2.435*initial_value while i&lt;(np.pi*(2.435*10**(27))): i=np.pi*2.435*initial_value[j] phi=[] phi_prime=[] for x in range (len(m)): def dzdt(z,T): return [z, -3*1.4441*(10**(-6))*m[x]*np.sqrt(0.69)*(mp.coth(1.5*np.sqrt(0.69)* (10**(-6))*1.4441*m[x]*T))*z - z] z0 = [i,0] ts = tp/m[x] zs = spi.odeint(dzdt, z0, ts) phi.append(zs[99,0]) phi_prime.append(zs[99,1]) phi_m[j]=phi phi_m_prime[j]=phi_prime j+=1 end_time=time.perf_counter() print(end_time-start_time,&quot;seconds&quot;) </code></pre> <p>I don't know what is the problem. how can I get correct results? or at least as accurate as possible? or maybe I should rewrite the code in another form? thank you.</p> <p>UPDATE_1: I increased the steps and decreased the step sizes in this way:</p> <pre><code>tp=[] t1=np.logspace(0,1,100) t2=np.logspace(1,3,500) t3=np.logspace(3,4,700) t4=np.logspace(4,5,800) t5=np.logspace(5,6,1000) </code></pre> <p>and used <code>tp.append()</code> for every single elements in t1,t2,...,t5 ; but there is still some false results for some &quot;m[x]&quot;'s like for x=8 to 13, surprisingly results for x&lt;8 and x&gt;13 are not too bad!</p> <p>UPDATE_2:</p> <p>once again I increased the steps number:</p> <pre><code>tp=[] t1=np.logspace(0,1,100) t2=np.logspace(1,3,5000) t3=np.logspace(3,4,7000) t4=np.logspace(4,5,8000) t5=np.logspace(5,6,10000) </code></pre> <p>and used <code>tp.append()</code> for the elements; so now I have &quot;tp&quot; with 20100 steps that has more steps and smaller step sizes in the parts that function has his rapid oscillations; after several hours it's still under running! I do not know if this method will help or not.</p> https://scicomp.stackexchange.com/q/41525 3 How is the surface Jacobian determinant calculated in FEM? Mohamed Abdelhamid https://scicomp.stackexchange.com/users/37417 2022-06-24T21:02:49Z 2022-06-30T16:12:44Z <p>I am currently trying to evaluate surface forces on a structure. I came across P356 in Bathe's Finite Element Procedures 2014 (example 5.8) in which he related the edge derivative from the global to natural coordinates through this equation:</p> <p><span class="math-container">$dl = det(\mathbf{J}^S) \ dr$</span></p> <p>where <span class="math-container">$r$</span> is the natural coordinate along the edge of interest. I am not sure how <span class="math-container">$\mathbf{J}^S$</span> or its determinant are calculated! I am not sure how to relate that to the edge length for example!</p> <p>Thanks in advance.</p> https://scicomp.stackexchange.com/q/41524 14 Why aren't Krylov subspace methods popular in the Machine Learning community compared to Gradient Descent? SPARSE https://scicomp.stackexchange.com/users/39167 2022-06-24T19:05:25Z 2022-06-30T11:35:42Z <p>Historically, iterative methods for solving relatively simple-structured systems <span class="math-container">$Ax=b$</span> with <span class="math-container">$A$</span> being a <span class="math-container">$4\times 4$</span> matrix or to find the eigenvalues of that matrix assuming in both problems that <span class="math-container">$A$</span> is non-singular were not quite popular and were not needed in pre-1950s as stationary methods or perhaps classical iterative algorithms such as Jacobi, Gauss-Seidel would handle the task easily.</p> <p>Then in the <span class="math-container">$1950$</span>s the rapid development of computational tools made it possible to deal with large and sparse matrix system (possibly a <span class="math-container">$100\times 100$</span> at that time in engineering applications). This is all thanks to the development of Krylov Subspace tools which at that time was considered to be &quot;high-tech&quot;. This is seen from various celebrated algorithms such as Arnoldi Iteration (Arnoldi-Lanczos 1950) and Conjugate Gradient (Stiefel,1952). During the 1970s and 1980s when we had tons and tons of transistors inside computers this made it possible for us to witness the birth of GMRES and BiCGSTAB which are known to handle relatively large and sparse systems.</p> <p>In fact, it was during that period that BLAS and LAPACK were developed and ever since have continued to upgrade various Krylov algorithms and achieve a high level of robustness.</p> <p>In machine learning, many problems such as regression, neural networks, and even recommender systems can one way or another be formulated as an optimization problem which in turn can be expressed as a linear system <span class="math-container">$Fx=z$</span>. The apparent default choice for solving such problems is the &quot;Gradient Descent Algorithm&quot; (and its derivatives) a popular optimizer algorithm that seems to be used literally everywhere in machine learning problems (movie recommender system, linear regression,...) even a famous variant is the ADAM Optimizer which in many deep learning libraries is set as the default optimizer.</p> <p>Given that numerical linear algebra libraries have been for the past 40 years enhancing Krylov algorithms, why haven't they made an impact in today's world of machine learning as they did in engineering communities?</p> <p>Let's say linear regression is the underlying machine learning problem, why is the gradient descent algorithm preferred over using BiCGSTAB with some restart when rewriting the regression problem as <span class="math-container">$Fx=z$</span> where <span class="math-container">$F$</span> could be sparse and <span class="math-container">$10000\times 10000$</span>?</p> <p>I believe one answer (though it may be biased) is that the majority of the machine learning community is not paying attention to the computational mechanism of the algorithm. For instance, you can go ahead and apply ADAM optimizer to solve a simple regression problem (overkill) or a complicated large regression problem and may end up witnessing your algorithm diverging without properly understanding why. I found a meme from medium.com that jokingly describes this:</p> <p><a href="https://i.stack.imgur.com/ilLzc.jpg" rel="noreferrer"><img src="https://i.stack.imgur.com/ilLzc.jpg" alt="enter image description here" /></a></p> <p><strong>Source:</strong> <a href="https://medium.com/analytics-vidhya/mechanism-of-gradient-descent-optimization-algorithms-7baed9c9c35e" rel="noreferrer">https://medium.com/analytics-vidhya/mechanism-of-gradient-descent-optimization-algorithms-7baed9c9c35e</a></p> <p><strong>Update 1</strong>: I would like to thank <code>davidhigh</code> for the detailed answer. I would like to add that my experience with Krylov methods were mainly focused on modelling dynamical system. Tim Davis has a <a href="https://sparse.tamu.edu/" rel="noreferrer">very interesting</a> collection of sparse matrices with detailed information on where the problem arises from (in particular engineering applications) Krylov methods are highly applicable there.</p> <p><strong>Update 2</strong>: <code>Neil Lindquist</code> made an important correction regarding BLAS and LAPACK not having Krylov solvers.</p> <p><strong>Final Remark</strong>: I welcome any additional answer especially if some users have experience with using Krylov methods in ML applications.</p> https://scicomp.stackexchange.com/q/41385 2 Place points at maximum distance in a convex 2D set Renan Andrade https://scicomp.stackexchange.com/users/43220 2022-05-23T13:32:10Z 2022-07-06T11:46:01Z <p>I need to place a given finite set of points within maximum distance of each other on 2D, constrained by a convex boundary, on Python.</p> <p>Honestly, I'm kind of lost. I have the explicit function to be minimized, but I'm not sure which algorithm I should use to adjust the points recursively. I thought of experimenting with clustering but I'm not sure if this would work.</p> <p>And, more importantly, I have no idea how to bound the space.</p> <p>Unfortunately, the convex boundary isn't a rectangle as well; it looks like this:</p> <p><a href="https://i.stack.imgur.com/N3ivq.png" rel="nofollow noreferrer"><img src="https://i.stack.imgur.com/N3ivq.png" alt="enter image description here" /></a></p> https://scicomp.stackexchange.com/q/40989 1 Global convergence behavior of several Krylov solvers in scipy.opt G. Fougeron https://scicomp.stackexchange.com/users/29718 2022-02-10T13:42:03Z 2022-07-03T15:35:35Z <p>In a context of mechanical simulation, where I solve the stationnary action principle directly (i.e. <span class="math-container">$\nabla S = 0$</span> for some scalar fonction <span class="math-container">$S$</span>), I use the wrapper scipy.optimize.newton_krylov to access several different non linear Krylov based solvers (in particular, &quot;lgmres&quot; and &quot;cgs&quot;).</p> <p>The equations have many solutions, and different initial guesses will return different solutions. This is to be expected. In order to find multiple solutions, I initialize the initial guess using a random number generator, and then call scipy.optimize.newton_krylov.</p> <p>My observation is the following: the two solvers &quot;lgmres&quot; and &quot;cgs&quot; tend to converge towards qualitatively different solutions. In a sense, the solutions found by &quot;cgs&quot; are nicer and simpler while those found by &quot;lgmres&quot; are more complex and convoluted.</p> <p>My question(s): Is there an explaination for this situation ? Is it possible that the dynamics of the &quot;cgs&quot; solver are unstable arround the most complex solutions, and why ? Do the basins of convergence of the two methods have qualitatively different properties ? And finally, is it possible to link this behavior to the spectrum of the Hessian of <span class="math-container">$S$</span> at the solution ?</p> <p>Thanks in advance,</p> https://scicomp.stackexchange.com/q/40120 0 How does one obtain tortoise coordinates by integrating with GNU Scientific Library (GSL)? JoséM https://scicomp.stackexchange.com/users/41359 2021-09-26T15:52:42Z 2022-07-04T07:03:25Z <p>I’m trying to to obtain the values of the tortoise coordinates (Eddington-Finkelstein Coordinates) integrating the expression: <span class="math-container">$\frac{dr^*}{dr} = (1 - rS/r)^{-1}$</span> using the GNU Scientific Library (GSL). I know that the analytical expression is <span class="math-container">$r* = r + rs·ln(r/rS - 1)$</span> so I can compare both results: analytical and integrated.</p> <p>Unfortunately, the integrated solution seems rotated</p> <p><a href="https://i.stack.imgur.com/wMLYQ.png" rel="nofollow noreferrer"><img src="https://i.stack.imgur.com/wMLYQ.png" alt="enter image description here" /></a></p> <p>but even if I rotate the graph, the result is shifted…</p> <p><a href="https://i.stack.imgur.com/eBESs.png" rel="nofollow noreferrer"><img src="https://i.stack.imgur.com/eBESs.png" alt="enter image description here" /></a></p> <p>If I add 2 constants the graphs are the same</p> <p><a href="https://i.stack.imgur.com/PqmLm.png" rel="nofollow noreferrer"><img src="https://i.stack.imgur.com/PqmLm.png" alt="enter image description here" /></a></p> <p>My question are:</p> <ol> <li>why the calculated is rotated?</li> <li>And why if I rotated the graph it’s still shifted? (with the values -x + 10.605, -yxCalculated[i] + rS) both graphs overlap)</li> </ol> <p>I need to know how to do this calculations because in the case of Black Hole Bounce or Wormholes we usually don't have an analytical expression.</p> <p>Here it’s my code:</p> <pre><code>#include &lt;stdio.h&gt; #include &lt;assert.h&gt; #include &lt;math.h&gt; #include &lt;gsl/gsl_errno.h&gt; #include &lt;gsl/gsl_odeiv2.h&gt; #define SYS_DIM 1 #define NODES 10001 #define SCH_MIN_INF -10000 #define SCH_PLU_INF 10000 struct yxParams { double rS; /* Schwarzschild radius */ }; double rS = 2; /* I refer to r* as y and r as x: r = x r* = y So yx = r*(r) and xy = r(r*) */ double yxAnalytical[NODES]; double yxCalculated[NODES]; /** @brief ODE IV 1st grade to solve y[x]' @param t Independent variable @param y[] Left side of the 1st grade system of equations @param sys[] Right side of the 1st grade system of equations */ int funcY(double t, const double y[], double sys[], void* params) { struct yxParams *par = (struct yxParams*)params; double rS = (par-&gt;rS); (void)(t); double x0 = y; assert(x0 != 0); sys = 1.0/(1.0 - rS/x0); return GSL_SUCCESS; } /** @brief r*(r) integration... so the inverse r(r*) tortoise coordinates. We do this because it's easier to get an analytic expression from r*(r) that from r(r*) @param IN, a: interval lower limit @param IN, b: interval upper limit @param IN, n: number of subintervals @param IN, ic: Initial condition @param OUT, yx[]: Array with the integration values */ int yx_integration(double a, double b, int n, double iC) { int status = GSL_SUCCESS; double h = (b - a)/n; double x0 = a; double x1 = x0 + h; double epsAbs = 0; double epsRel = 1e-6; struct yxParams params = {rS}; gsl_odeiv2_system sys = {funcY, NULL, SYS_DIM, &amp;params}; const gsl_odeiv2_step_type* T = gsl_odeiv2_step_rkf45; gsl_odeiv2_driver* d = gsl_odeiv2_driver_alloc_y_new (&amp;sys, T, h, epsAbs, epsRel); double y = {iC}; int i; for (i = 0; i &lt; NODES; i++) { status = gsl_odeiv2_driver_apply(d, &amp;x0, x1, y); x0 = x1; x1 = x0 + h; if (status != GSL_SUCCESS) { printf (&quot;error, return value = %d\n&quot;, status); break; } else { double tmp = y; if(_isnanf(tmp) || isinf(tmp)) { tmp = SCH_MIN_INF; } yxCalculated[i] = tmp; } } gsl_odeiv2_driver_free(d); return status; } /** @brief Tortoise analytical expression: r*(r) @param IN a: interval lower limit @param IN h: step */ void yx_analytical(double a, double h) { double x = a; int i; for(i = 0; i &lt; NODES; i++) { double tortoise = x + rS*log(x/rS - 1); if(_isnanf(tortoise)) { tortoise = SCH_MIN_INF; } yxAnalytical[i] = tortoise; x += h; } } /** @brief Save the data so we can plot it @param IN x @param IN h */ void txt_data(double x, double h) { FILE *fp; fp = fopen(&quot;data.txt&quot;, &quot;w&quot;); fprintf(fp, &quot;#%32s %32s %32s %32s\n&quot;, &quot;yxAnalytical&quot;, &quot;x&quot;, &quot;yxCalculated&quot;, &quot;x&quot;); int i; for (i = 0; i &lt; NODES; i++) { /* If we plot 'yx vs x', we're plotting the inverse function... so we're plotting r(r*) the tortoise coordinates xy[]. Q1) yxCalculated is rotated with respect yxAnalytical. Why? Q2) if I plot (if I'm not wrong this implies rotating 180º the graph) fprintf(fp, &quot;%.32f %.32f %.32f %.32f\n&quot;, yxAnalytical[i], x, -x, -yxCalculated[i]); both graphs are practically the same... but not equal. I need to shift the graph fprintf(fp, &quot;%.32f %.32f %.32f %.32f\n&quot;, yxAnalytical[i], x, -x + 10.605, -yxCalculated[i] + rS); so calculations seems to be in the right direction... but certainly something is really wrong */ fprintf(fp, &quot;%.32f %.32f %.32f %.32f\n&quot;, yxAnalytical[i], x, yxCalculated[i], x); /* Trick to overlap the graphs... * fprintf(fp, &quot;%.32f %.32f %.32f %.32f\n&quot;, yxAnalytical[i], x, -x + 10.605, -yxCalculated[i] + rS); */ x += h; } fclose(fp); } /**----------------------------------------------------------------------------- ------------------------------------------------------------------------------*/ int main() { /* To our purposes these values are enough far away from the BH. These values are for the numerical integration. */ double a = -150; double b = 150; int n = NODES - 1; double h = (b - a)/n; /* Do calculations... */ yx_analytical(a, h); yx_integration(a, b, n, a); /* Save the data in a text file */ txt_data(a, h); return 0; } </code></pre> <p><strong>I'll answer myself</strong></p> <p>This's my final code (integrating backwards). I've added the jacobian and solve a little issue (misconception) with my funcY.</p> <pre><code>#include &lt;stdio.h&gt; #include &lt;assert.h&gt; #include &lt;math.h&gt; #include &lt;gsl/gsl_errno.h&gt; #include &lt;gsl/gsl_odeiv2.h&gt; #include &lt;gsl/gsl_matrix.h&gt; #define SYS_DIM 1 #define NODES 10001 #define SCH_MIN_INF -10000 #define SCH_PLU_INF 10000 struct yxParams { double rS; /* Schwarzschild radius */ }; double rS = 2; /* I refer to r* as y and r as x: r = x r* = y So yx = r*(r) and xy = r(r*) */ double yxAnalytical[NODES]; double yxCalculated[NODES]; /** @brief yx[] jacobian, so you can use implicit methods rS/(r^2*(1 - rS/r)^2) */ int jacoY(double r, const double y[], double *dfdy, double dfdt[], void *params) { (void)(y); struct yxParams *par = (struct yxParams*)params; double rS = (par-&gt;rS); gsl_matrix_view dfdy_mat = gsl_matrix_view_array(dfdy, 1, 1); gsl_matrix *m = &amp;dfdy_mat.matrix; /*gsl_matrix_set (m, 0, 0, -rS/(r*r*(rS/r - 1)*(rS/r - 1)));*/ gsl_matrix_set (m, 0, 0, -rS/(r*r*(1 - rS/r)*(1 - rS/r))); dfdt = 0.0; return GSL_SUCCESS; } /** @brief ODE IV 1st grade to solve y[x]' @param x Independent variable @param y[] Left side of the 1st grade system of equations @param sys[] Right side of the 1st grade system of equations */ int funcY(double x, const double y[], double sys[], void* params) { struct yxParams *par = (struct yxParams*)params; double rS = (par-&gt;rS); assert(x != 0); sys = 1.0/(1.0 - rS/x); return GSL_SUCCESS; } /** @brief r*(r) integration... so the inverse r(r*) tortoise coordinates. We do this because it's easier to get an analytic expression from r*(r) that from r(r*) @param IN, a: interval lower limit @param IN, b: interval upper limit @param IN, n: number of subintervals @param IN, ic: Initial condition @param OUT, yx[]: Array with the integration values */ int yx_integration(double a, double b, int n, double iC) { int status = GSL_SUCCESS; double h = (b - a)/n; double x0 = a; double x1 = x0 + h; double epsAbs = 0; double epsRel = 1e-6; struct yxParams params = {rS}; gsl_odeiv2_system sys = {funcY, jacoY, SYS_DIM, &amp;params}; const gsl_odeiv2_step_type* T = gsl_odeiv2_step_rk8pd; /* rkf45, rkck, rk8pd; rk4imp bsimp msadams msbdf */ gsl_odeiv2_driver* d = gsl_odeiv2_driver_alloc_y_new (&amp;sys, T, h, epsAbs, epsRel); double y = {iC}; int i; for (i = 0; i &lt; NODES; i++) { status = gsl_odeiv2_driver_apply(d, &amp;x0, x1, y); x0 = x1; x1 = x0 + h; if (status != GSL_SUCCESS) { printf (&quot;error, return value = %d\n&quot;, status); printf(&quot;x0: %f, y: %f&quot;, x0, y); break; } else yxCalculated[i] = y; } gsl_odeiv2_driver_free(d); return status; } /** @brief Tortoise analytical expression: r*(r) @param IN a: interval lower limit @param IN h: step */ void yx_analytical(double a, double h) { double x = a; int i; for(i = 0; i &lt; NODES; i++) { double tortoise = x + rS*log(x/rS - 1); if(isnan(tortoise)) tortoise = SCH_MIN_INF; yxAnalytical[i] = tortoise; x += h; } } /** @brief Save the data so we can plot it @param IN x @param IN h */ void txt_data(double x, double h) { FILE *fp; fp = fopen(&quot;data.txt&quot;, &quot;w&quot;); int n = NODES - 1; int i; for (i = 0; i &lt; NODES; i++) { fprintf(fp, &quot;%.32f %.32f %.32f %.32f\n&quot;, yxAnalytical[i], x, yxCalculated[n - i], x); x += h; } fclose(fp); } /**----------------------------------------------------------------------------- ------------------------------------------------------------------------------*/ int main() { /* To our purposes these values are enough far away from the BH. These values are for the numerical integration. */ double a = -150; double b = 150; int n = NODES - 1; double h = (b - a)/n; /* Do calculations... */ yx_analytical(a, h); /*yx_integration(b, a, n, b);*/ /* Unfortunately we need and exact IC... we cannot use the tortoise approximation rS = r enough far away... */ double ic = b + rS*log(b/rS - 1); yx_integration(b, a, n, ic); /* Save the data in a text file */ txt_data(a, h); return 0; } </code></pre> <p>So if we use the tortoise approximation for the initial condition, where rS = r far away, we've got a curve with a perfect shape... but with an offset <a href="https://i.stack.imgur.com/riymT.png" rel="nofollow noreferrer"><img src="https://i.stack.imgur.com/riymT.png" alt="enter image description here" /></a> note that at the horizon we cannot continue with the integration (but is ok). If we use the exact IC, where <span class="math-container">$r_* = r + rS·ln(r/rS - 1)$</span> we get the exact curve <a href="https://i.stack.imgur.com/7KY2W.png" rel="nofollow noreferrer"><img src="https://i.stack.imgur.com/7KY2W.png" alt="enter image description here" /></a> except that again, we cannot continue the integration at the horizon.</p> https://scicomp.stackexchange.com/q/35600 5 4th order tensor rotation - sources to refer Sagar Shah https://scicomp.stackexchange.com/users/36594 2020-07-23T20:43:34Z 2022-07-06T01:27:34Z <p>I am trying to model a linear elastic material in Abaqus using a UMAT. For my application, I need to rotate the 6x6 compliance matrix for a given set of eigenvectors (or a rotation matrix). I came across a thread titled <a href="https://scicomp.stackexchange.com/questions/8093/debugging-a-rotation-matrix-for-elastic-constants">&quot;debugging a rotation matrix for elastic constants&quot;</a> where this theory was explained in very good detail.</p> <p>I was wondering if there are any good sources out there that show the actual matrix rotation process either in a computational sense or just a theoretical representation.</p> https://scicomp.stackexchange.com/q/20644 10 Use of machine learning in computational fluid dynamics EngrStudent https://scicomp.stackexchange.com/users/3979 2015-09-04T11:22:27Z 2022-07-01T12:28:36Z <p><strong>Background:</strong><br /> I have only built one working numeric solution to 2d Navier-Stokes, for a course. It was a solution for lid-driven cavity flow. The course, however, discussed a number of schemas for spatial discretizations and time discretizations. I have also taken more symbol-manipulation coursework applied to NS.</p> <p>Some of the numeric approaches to handle conversion of the analytic/symbolic equation from PDE to finite difference include:</p> <ul> <li>Euler FTFS, FTCS, BTCS</li> <li>Lax</li> <li>Midpoint Leapfrog</li> <li>Lax-Wendroff</li> <li>MacCormack</li> <li>offset grid (spatial diffusion allows information to spread)</li> <li>TVD</li> </ul> <p>To me, at the time, these seemed like &quot;insert-name finds a scheme and it happens to work&quot;. Many of these were from before the time of &quot;plentiful silicon&quot;. They are all approximations. In the limit they. in theory, lead to the PDE's.</p> <p>While Direct Numerical Simulation (<a href="https://en.wikipedia.org/wiki/Direct_numerical_simulation" rel="nofollow noreferrer">DNS</a>) is fun, and Reynolds Averaged Navier-Stokes (<a href="https://en.wikipedia.org/wiki/Reynolds-averaged_Navier%E2%80%93Stokes_equations" rel="nofollow noreferrer">RANS</a>) is also fun, they are the two &quot;endpoints&quot; of the continuum between computationally tractable, and fully representing the phenomena. There are multiple families of approaches that live interior to these.</p> <p>I have had CFD professors say, in lecture, that most CFD solvers make pretty pictures, but for the most part, those pictures do not represent reality and that it can be very tough, and take lots of work, to get a solver solution that does represent reality.</p> <p>The sequence of development (as I understand it, not exhaustive) is:</p> <ol> <li><p>start with the governing equations -&gt; PDE's</p> </li> <li><p>determine your spatial and temporal discretization -&gt; grid and FD rules</p> </li> <li><p>apply to the domain including initial conditions and boundary conditions</p> </li> <li><p>solve (lots of variations on matrix inversion)</p> </li> <li><p>perform gross reality checks, fit to known solutions, etc..</p> </li> <li><p>build some simpler physical models derived from analytic results</p> </li> <li><p>test them, analyze, and evaluate</p> </li> <li><p>iterate (jumping back to either step 6, 3, or 2)</p> </li> </ol> <p>Thoughts:<br /> I have recently been working with CART models, oblique trees, random forests, and gradient boosted trees. They follow more mathematically derived rules, and the math drives the shape of the tree. They work to make discretized forms well.</p> <p>Although these human-created numeric approaches work somewhat, there is extensive &quot;voodoo&quot; needed to connect their results to the physical phenomena they are meant to model. Often the simulation does not substantially replace real-world testing and verification. It is easy to use the wrong parameter, or not account for variation in geometry or application parameters experienced in the real world.</p> <p><strong>Questions:</strong></p> <ul> <li>Has there been any approach to let the nature of the problem define<br /> the appropriate discretization, spatial and temporal differencing scheme, initial conditions, or solution?</li> <li>Can a high definition solution coupled with the techniques of machine learning be used to make a differencing scheme that has much larger step sizes but retains convergence, accuracy, and such?</li> <li>All of these schemes are accessibly &quot;humanly tractable to derive&quot; - they have a handful of elements. Is there a differencing scheme with thousands of elements that does a better job? How is it derived?</li> </ul> <p>Note: I will follow up with the empirically intialized and empirically derived (as opposed to analytically) in a separate question.</p> <p><strong>UPDATE:</strong></p> <ol> <li><p>Use of deep learning to accelerate lattice Boltzmann flows. Gave ~9x speedup for their particular case</p> <p>Hennigh, O. (in press). Lat-Net: Compressed Lattice Boltzmann Flow Simulations using Deep Neural Networks. Retrieved from: <a href="https://arxiv.org/pdf/1705.09036.pdf" rel="nofollow noreferrer">https://arxiv.org/pdf/1705.09036.pdf</a></p> <p>Repo with code (I think):<br /> <a href="https://github.com/loliverhennigh/Phy-Net" rel="nofollow noreferrer">https://github.com/loliverhennigh/Phy-Net</a></p> </li> <li><p>About 2 orders of magnitude faster than GPU, 4 orders of magnitude, or ~O(10,000x) faster than CPU, and same hardware.</p> <p>Guo, X., Li, W. &amp; Ioiro, F. Convolutional Neural Networks for Steady Flow Approximation. Retrieved from: <a href="https://autodeskresearch.com/publications/convolutional-neural-networks-steady-flow-approximation" rel="nofollow noreferrer">https://autodeskresearch.com/publications/convolutional-neural-networks-steady-flow-approximation</a></p> </li> <li><p>Others who have looked into the topic about 20 years ago:</p> <p>Muller, S., Milano, M. &amp; Koumoutsakos P. Application of machine learning algorithms to flow modeling and optimization. Center for Turbulence Research Annual Research Briefs 1999 Retrieved from: <a href="https://web.stanford.edu/group/ctr/ResBriefs99/petros.pdf" rel="nofollow noreferrer">https://web.stanford.edu/group/ctr/ResBriefs99/petros.pdf</a></p> </li> </ol> <p><strong>Update (2017):</strong><br /> <a href="https://arxiv.org/abs/1712.06567" rel="nofollow noreferrer">This</a> characterises the use of non-gradient methods in deep learning, an arena which has been exclusively gradient based. While the direct implication of activity is in deep learning, it also suggests that GA can be used as an equivalent in solving a very hard, very deep, very complex problem at the level consistent with or superior to gradient descent based methods.</p> <p>Within the scope of this question, it might suggest that a larger-scale, machine-learning based attack might allow &quot;templates&quot; in time and space that substantially accelerate convergence of gradient-domain methods. The article goes as far as to say that sometimes going in the direction of gradient descent moves away from the solution. While in any problem with local optima or pathological trajectories (most high-value real-world problems have some of these) it is expected that the gradient isn't globally informative, it is still nice to have it quantified and validated empirically as it was in this paper and the ability to &quot;jump the bound&quot; without requiring &quot;reduction of learning&quot; as you get in momentum or under-relaxation.</p> <p><strong>Update (2019):</strong><br /> It seems that google now has a contribution &quot;how to find a better solver&quot; piece of the AI puzzle. <a href="https://ai.googleblog.com/2018/08/introducing-new-framework-for-flexible.html" rel="nofollow noreferrer">link</a> This is a part of making the AI make the solver.</p> <p><strong>Update (2020):</strong><br /> And now they are doing it, and doing it well...<br /> <a href="https://arxiv.org/pdf/1911.08655.pdf" rel="nofollow noreferrer">https://arxiv.org/pdf/1911.08655.pdf</a></p> <p>Is that a 100x speedup for Navier-Stokes solution??</p> <p>It could be argued that they could then deconstruct their NN to determine the actual discretization. I particularly like figure 4.</p> <p>Update 2020 <a href="https://arxiv.org/pdf/2006.08762.pdf" rel="nofollow noreferrer">link</a></p> <p><strong>Update 2022</strong> - Physics Informed neural networks allow state compression by a factor of 85x with only slight loss. (<a href="https://arxiv.org/pdf/2201.03617v2.pdf" rel="nofollow noreferrer">link</a>) One of the values of being able to stop, store the state, move it to another computer, and be able to pick it up and restart at the stop-point and not the initial condition is that it can allow reproducibility. A decent near-estimate of a state means the compute engine can get to the actual with only a little work. A ML-based &quot;picture&quot; that is good enough can act just like all the compute overhead from the initial state to the good enough estimate.</p> https://scicomp.stackexchange.com/q/8125 11 Are there any libraries out there that implement block Krylov subspace methods? Geoff Oxberry https://scicomp.stackexchange.com/users/276 2013-08-01T01:40:19Z 2022-07-04T03:37:44Z <h1>Question</h1> <p>Are there libraries out there that implement block Krylov subspace methods? (I was not able to find any from a simple Google search.)</p> <h1>Background</h1> <p>Right now, I am working with a code that solves several systems of the form</p> <p>\begin{align} Ax_{i} = b_{i}, \end{align}</p> <p>where $A$ is $n$ by $n$, $i = 1, \ldots, m$, and in general, $n \gg m$ ($n$ is at least 10 times larger than $m$, and could be several orders of magnitude larger than $m$). In addition, $A$ is large enough that iterative methods are preferable.</p> <p>Currently, the code solves $Ax_{i} = b_{i}$ with a Krylov subspace method for each separate right-hand side with a given preconditioner, and this approach seems to work well enough. (The code also recomputes the preconditioner, which is unnecessary.)</p> <p>To speed up the code, it seems like it could be worth looking at solving</p> <p>\begin{align} AX = B \end{align}</p> <p>instead, where $B = [b_{1}\, b_{2}\, \ldots\, b_{m}]$ is an $n$ by $m$ matrix gathering up all of the $b_{i}$, and $X$ is also an $n$ by $m$ matrix gathering up the $x_{i}$.</p> <p>I'd like to try out solving $AX = B$ with block Krylov subspace methods (using the same preconditioner as before) to see if it outperforms the current approach of solving $Ax_{i} = b_{i}$ and looping over $i$. It is both possible and unlikely that the solution to $Ax_{i} = b_{i}$ could then be used as a guess for the system $Ax_{j} = b_{j}$ for $i \neq j$; I don't expect there to be a relationship among the various systems.</p> https://scicomp.stackexchange.com/q/2630 12 Libraries for solving Lyapunov's equation NRH https://scicomp.stackexchange.com/users/44 2012-06-26T15:28:16Z 2022-06-30T19:19:38Z <p>The following matrix equation $$B\Sigma + \Sigma B^T + C = 0$$ in $\Sigma$ $-$ for given $B$ and $C$ matrices $-$ appears in my work as a characterization of a covariance matrix. I have learned that this equation is known, in particularly in continuous time control theory, as <strong>Lyapunov's equation</strong>, and that there are various well known algorithms for solving it that exploit the special nature of this linear equation. </p> <p>From googling I have also learned that there exist Matlab and Fortran implementations. I have found SLICOT and RECSY. Due to licensing issues access to SLICOT source has been stopped, though. </p> <p>Most of my work is implemented in R, and as I have been unable to find an R interface to a solver, I consider writing one myself. My question is then if SLICOT is the best available Fortran (or C) library with an implementation of a solver of Lyapunov's equation? I am also interested in implementations that can handle large sparse $B$ matrices. </p> https://scicomp.stackexchange.com/q/2490 4 Line search for Newton method Alexander https://scicomp.stackexchange.com/users/1089 2012-06-11T10:29:01Z 2022-06-29T16:00:27Z <p>If we want to solve nonlinear minimization problem </p> <p>$$\min_{x} f(x),$$</p> <p>making least-squares assumption and using Gauss-Newton method so that at k$th$ iteration we have:</p> <p>$$J_k^T J_k p_k = - J_k^T r_k,$$</p> <p>with vector $r$ being residual and matrix $J$ Jacobian.</p> <p>We then update $x$ in the following way:</p> <p>$$x_{k+1} = x_k + \alpha p_k,$$</p> <p>where $\alpha \in (0, 1]$</p> <p>The question is how to find $\alpha$ such that $f(x_{k+1}) = min$</p> <p>In other words, what step length calculating strategy would be good enough taking into account $f$ and $f'$ are expensive to compute and highly nonlinear?</p> <p>I'm aware of methods that approximate this problem with polynomial, e.g. in case of quadratic approximation:</p> <p>$$p_0 + p_1 \alpha + p_2 \alpha^2 = min$$</p> <p>where $p_0 = f(x_k), p_1 = f'(x_k), p_2 = f(x_k + \alpha p_k)$</p> <p>But I'm wondering what are the other options to try? Can somebody point me to a good overview or shortly write down different techniques. </p>