numerical-analysis

[fresh oxide]

ok cool

faster than row reduce even?

i dont quite understand the proof tbh

[mint hemlock]

I’m not 100% sure what is faster in this scenario, I’m opted to think power method still

O(n^3) is just icky, but so is O(n^2*k)

[fresh oxide]

so essentially power method converges faster when?

the 1 - 2nd largest eigenvalue is high?

@mint hemlock i assume this is still faster?

[mint hemlock]

Oh yes, that’s one of the methods I was forgetting

There may be a way to leverage knowing that 1 is an eigenvalue

[fresh oxide]

thats what i was thinking

btw dot product of two matrices is o(n^3) right? @mint hemlock

[mint hemlock]

yes

There are interesting algorithms that make it less, but the O(n^3) approach is basically what you get

[fresh oxide]

im thinking of trying to justify using power method 😂

[mint hemlock]

why would you be doing matrix products there

[fresh oxide]

nah that was for a different thing

just wondering

[mint hemlock]

ah okay

[fresh oxide]

its RNN ANN LSTM models

def dot(x, y):
    assert x.shape[1] == y.shape[0], "Incorrect matrix dimensions for dot product. The number of columns of the 1st matrix must equal to the number of rows of the 2nd matrix."
    dotProduct = np.zeros((x.shape[0], y.shape[1]))
    for i in range(0, x.shape[0]):
        for j in range(0, y.shape[1]):
            for k in range(0, x.shape[1]):
                dotProduct[i][j] += x[i][k] * y[k][j]
    return dotProduct


def transpose(x):
    transposedMat = np.zeros((x.shape[1], x.shape[0]))
    for i in range(0, x.shape[0]):
        for j in range(0, x.shape[1]):
            transposedMat[j][i] = x[i][j]
    return transposedMat


def hadamard(x, y):
    assert x.shape[1] == y.shape[1] and x.shape[0] == y.shape[0], "Incorrect matrix dimensions for hadamard product. The dimensions of both matrices must be equal for an elementwise product."
    hadamardProduct = np.zeros((x.shape[1], x.shape[0]))
    for i in range(0, x.shape[0]):
        for j in range(0, x.shape[1]):
            hadamardProduct[i][j] = x[i][j] * y[i][j]
    return hadamardProduct


def outer(x, y):
    return dot(transpose(x), y)

i wanna know if this code is correct

[mint hemlock]

it appears correct

I can’t test it against anything because I’m not at a computer

[fresh oxide]

hm ok

so matrix dot product is indeed o(n^3)?

[mint hemlock]

okay so multiplication operation is o(n^3)
You can perform the dot product operation without doing transpose(x)

[prime kraken]

for square matrices, the pen and paper method is indeed O(n^3)

that outer product definition is only for vectors tho

everything looks ok

[fresh oxide]

wdym?

[fresh oxide]

ok thnx 👍

i dont quite get this proof tbh:

Here's a really elementary proof (which is a slight modification of Fanfan's answer to a question of mine). As Calle shows, it is easy to see that the eigenvalue 1 is obtained. Now, suppose Ax=λx for some λ>1. Since the rows of A are nonnegative and sum to 1, each element of vector Ax is a convex combination of the components of x, which can be no greater than xmax, the largest component of x. On the other hand, at least one element of λx is greater than xmax, which proves that λ>1 is impossible.

Theorem: The largest eigenvalue of a stochastic matrix is 1.

Proof: First, if A is a stochastic matrix, then A1 = 1, since each row of A sums to 1. This proves that 1 is an eigenvalue of A. Second, suppose there exists λ > 1 and nonzero x such that Ax = λx. Let xi be a largest element of x. Since any scalar multiple of x will also satisfy this equation we can assume, without loss of generality, that xi > 0. Since the rows of A are nonnegative and sum to 1, each entry in λx is a convex combination of the elements of x. Thus no entry in λx can be larger than xi. But since λ > 1, λxi > xi. Contradiction. Therefore, the largest eigenvalue of A is 1.

[mint hemlock]

transpose(A) isn’t necessary for matrix multiplication is all I meant

[fresh oxide]

u mean in outer or in dot?

if in dot, im not sure how i would code it otherwise?

[mint hemlock]

in either case, the transpose of a matrix just flips [a][b] to [b][a], so you can just multiply along that instead. The transpose operation in your code is O(n^2)

This is just a minor overall point that is just implementation related

you’re ideas are correct

[fresh oxide]

so u mean i can do dot product in O(n^2)?

instead of O(n^3) like it is rn?

[mint hemlock]

no

Your implementation is fine

It’s just not as efficient as it could be is all

O(n^3(multiplication)+n^2(transpose)) is still O(n^3)

again, this is just more of a nitpick on implementation that i’m just referring to

[fresh oxide]

u mean for outer product right?

not for dot product?

[mint hemlock]

Yes I guess

[fresh oxide]

ok thnx

[mint hemlock]

the idea in this proof is that because each row adds to 1,
if $\lambda > 1$, it would be lead to contradiction because $\lambda x_m = a_{m,1}x_1+…+a_{m,n}x_n > x_m$ for every m
in other words, $\lambda x_m \le max(x_i)$ [max element of x]
but $\lambda max(x_i) > 1$, which is a contradiction

[fresh oxide]

ah ok

so sum will be greater than 1 basically?

[pine jetty]

kirby

[mint hemlock]

yeah, so if the largest eigenvalue were to be more than 1, there would be an issue because the rows of A each add up to 1

and of course, the upper bound of lambda * x_m for each m is the maximum element of x (because the rows of A add up to 1)

[fresh oxide]

ah ok thnx

[fresh oxide]

i dont quite get this tbh

[mint hemlock]

so $x =[x_1,x_2,…,x_n]^T$ such that $Ax = \lambda x,, \lambda > 1$
wlog let$ x_1 > x_2 > … > x_n$
then [(Ax)1 = A{1,1}x_1 + A_{2,1}x_2 + … + A_{n,1}x_n]
However, every $x_2,x_3,…,x_n < x_1$ so
[(Ax)1 \le x_1(A{1,1}+A_{2,1} +…+ A_{n,1})]
But $A_{1,1}+A_{2,1}+\dots+A_{n,1} =1$ so $(Ax)_1\le x_1$, but $(Ax)_1 = \lambda x_1 > x_1$, so it’s a contradiction.

[pine jetty]

kirby

[mint hemlock]

\lambda is the maximal eigenvalue

[onyx snow]

I noticed when using the numpy.fft.fft and numpy.fft.fftfreq there always seems to be a peak at x=0, even if that is not a sample frequency. I'm wondering if this is due to the fourier transform functions or something mathematical with the fourier transform itself?

[mint hemlock]

@onyx snow All I can think of is checking if there is a negative offset or something in the signal you're using. A potential solution I found is to subtract the average value from the signal in the time domain

Yeah, it looks like the frequency at 0Hz is just the mean of the data

so before performing fft, subtract the average of the data and see if there's still an issue @onyx snow

[onyx snow]

Ah perfect thank you!

[mint hemlock]

Yeah, it doesn't make too much sense in a mathematical sense, but indeed would make sense in context to there being some sort of bias

[wide spear]

I think it makes a lot of mathematical sense though

The 0 frequency component is the mean and everything else oscillates around it

You get the same thing with Fourier transforms

[wide spear]

Convolution with a smooth function

You convolute f with a Gaussian

[prime kraken]

you could also look at the expressions for butterworth filters, which approximate this type of function

they're of the form log(ratio of polynomials)

otherwise, what angetenar said will work, but you'd possibly want to use a 0-phase filter so that you don't get a sort extra flat part at the beginning

[prime kraken]

that's essentially the same thing

you want the derivative to be zero in the interval?

[wide spear]

You fill in f with whatever you want from 1 to 2

[prime kraken]

i don't get what you mean by this

[wide spear]

(texit is currently broken)

[prime kraken]

use splines of the order you like

why not?

splines are defined so that the derivatives match at the endpoints

yes, splines are defined by making a system of equations out of the function values and derivatives at the endpoints

it doesn't matter what happens in the interval (1,2)

the derivatives at 1 and 2 will match those of the other components of the piecewise function

probably needs to be at least cubic

maybe order 5 and use the endpoints and first 2 derivatives

that'll give you 6 equations

do you need it to be really infinitely differentiable?

and do you need the expression to be defined analytically

i can't think of any way off the top of my head, and i don't think it can be done with polys

[pine jetty]

squirtlespoof

squirtlespoof

[edgy sage]

@lilac forge Is $f_n$ short for $f(t_n, y_n)$?

[pine jetty]

IlIIllIIIlllIIIIllll

[edgy sage]

@lilac forge Find the largest $n$ for which the BDF2 $f(0) \approx a_1f(-2h) + a_2f(-h) + a_3hf'(0)$ is exact on $f(x) = x^n$.

[pine jetty]

IlIIllIIIlllIIIIllll

[edgy sage]

It's just a linear system, and I think it gives you the coefficients you wrote

@lilac forge No, to derive the coefficients

Yeah but once you know that the approximation is exact on $1, x, x^2$, you get the LTE for free

[pine jetty]

IlIIllIIIlllIIIIllll

[edgy sage]

Through Taylor's theorem

Let $E(f) = f(0) - (a_1f(-2h) + a_2f(-h) + a_3hf'(0))$. Note $E$ is linear. We have $$f(x) = f(0) + f'(0)x + \frac{f''(0)}{2}x^2 + \frac{f^{(3)}(0)}{3!}x^3 + \frac{f^{(4)}(\theta(x)x)}{4!}x^4$$

[pine jetty]

IlIIllIIIlllIIIIllll

[edgy sage]

Take $E$ to get $E(f) = \frac{f^{(3)}(0)}{3!}E(x^3) + E(\frac{f^{(4)}(\theta(x)x)}{4!}x^4)$.

[pine jetty]

IlIIllIIIlllIIIIllll

[edgy sage]

@lilac forge It's because that's how the coefficents were found

@lilac forge You find the coefficients by mandating that the approximation be exact on $1, x, x^2$. This gives you a linear system, which has a unique solution

[pine jetty]

IlIIllIIIlllIIIIllll

squirtlespoof

[edgy sage]

I would do $$f(x) = f(0) + f'(0)x + \frac{f''(0)}{2}x^2 + \frac{f^{(3)}(0)}{3!}x^3 + \frac{f^{(4)}(\theta(x)x)}{4!}x^4$$

[pine jetty]

IlIIllIIIlllIIIIllll

[edgy sage]

@dire wren I forgot to say I am translating coordinates so that $t_n = 0$

[pine jetty]

IlIIllIIIlllIIIIllll

[edgy sage]

and my $f$ is not the $f$ in the ode

[pine jetty]

IlIIllIIIlllIIIIllll

[edgy sage]

yeah

Now apply $E$ on both sides

[pine jetty]

IlIIllIIIlllIIIIllll

[edgy sage]

You get $E(f) = \frac{f^{(3)}(0)}{3!}E(x^3) + E(\frac{f^{(4)}(\theta(x)x)}{4!}x^4)$.

[pine jetty]

IlIIllIIIlllIIIIllll

[edgy sage]

$E(x^3)$ is simple to compute. As long as $f^{(4)}$ is bounded, the $E(\frac{f^{(4)}(\theta(x)x)}{4!}x^4)$ is $O(h^4)$, so it can be ignored if you are only interested in the leading error term.

[pine jetty]

IlIIllIIIlllIIIIllll

[split rain]

Hey I was wondering if anyone could help me with this problem:

I have that $r = -Ex$, but I can only seem to show that $||r||_2 \leq ||E||_F$

[pine jetty]

skwol

[edgy sage]

@split rain Let $E$ be defined by $Ex = -r$ and $E = 0$ on $\text{span}(x)^{\perp}$.

[pine jetty]

IlIIllIIIlllIIIIllll

[rich hatch]

Anyone have good resources on how to do multivariate pade approximants? From my googling around, there seem to be several methods that all have different strengths and weaknesses, and I feel like I'm lost in a bit of a twisty maze. Is there a location with an easy run down of my options?

[pine jetty]

TheRedLotus

[wide spear]

Code questions go in #computing-software

[uneven agate]

Moved it over. Sorry for the confusion.

[tall solar]

Anybody have any insight as to why both the SVD and the DFT can be used for compression?

It's like in both scenarios the flowchart it's

1. basis transformation
2. Thresholding
3. Inverse basis transformation

That's so weird to me. I've heard parsevals theorem used to explain why the dft can compress data but I'm still unsure about how the SVD fits into this

[prime kraken]

in the DFT case, you represent an image in a basis of sinusoids and throw away the ones that have small amplitudes, which is often the high frequency ones for natural images

the SVD one does not have such a nice physical interpretation, but is optimal in some sense

the SVD decomposes a matrix into a sum of rank 1 matrices that are made by multiplying the outer product of 2 unit norm vectors by a scaling factor

one then throws away the rank 1 matrices whose scaling factor is below some threshold

this is optimal in the sense that the result is the best rank K approximation to the original matrix in the Frobenius norm sense

this is the eckart-young-mirsky theorem

[tall solar]

Hmmm thanks for that explaination. I wonder if there is a similar optimality result for the dft?

Something like Fourier rank k "this is the most optimal function who's Fourier transform is a k-sparse vector"

[prime kraken]

i'm not familiar with results to that effect. not to say they don't exist, just really that i haven't run into them, for whatever reason

[tall solar]

Thank you 😁

[prime kraken]

well

in a more general direction

if the representation of the image is exactly K sparse in a DFT basis, then that's that

but that is usually not the case for images

in more general signal processing, sums of harmonics tend to show up as plausible models and then one does do this

then it becomes a question of a choice of basis for which the K sparsity yields the minimum error of some kind

[fathom rain]

And in some cases they are the same thing (e.g., circulant matrices)

[prime kraken]

oh that's a good example where the two match.

[tall solar]

What can I Google to learn more about this?

[prime kraken]

in signal processing/estimation theory this goes in the direction of "dictionary learning"

where you possibly have a known level of sparsity and some norm you want to minimize, and the goal is to learn the matrix, which may represent either a basis or an overcomplete set of vectors with which some other vector is computed

[tall solar]

Hmm now that's getting closer to the language I'm trying to speak. Yeah it has something to do with dictionary learning and sparse representations. I wonder if there is a general theory.

[prime kraken]

i'm pretty sure you find tons of papers on google scholar looking up dictionary learning

it's also closely related to blind learning/blind deconvolution

[near trench]

Hey guys, I have a question that says prove that a language $L$ is decidable iff $L \leq_m L(000^(11+111))$. So I think it's basically saying that prove a language $L$ is decidable iff it reduces to the language decided by the deterministic finite automaton based on the regular expression $000^(11+111)$

[pine jetty]

darkninja175

[near trench]

How does that work? One direction (<=) of implication is super easy as I can point out/prove that L(000*(11+111)) is decidable, and since L reduces to it then it must also be decidable

but the other implication says that any general decidable language can be reduced to L(000*(11+111)), which doesn't make sense to me at all

This is a computability/computer language course.

[edgy sage]

@near trench what is $\leq_m$?

[pine jetty]

IlIIllIIIlllIIIIllll

[near trench]

A language reduction

https://en.wikipedia.org/wiki/Turing_reduction

intuitively it's like reducing a problem down to a simpler equivalent problem. In this case it's reducing general decidability of a language. Formally it implies the existence of a reduction function that carries certain structure-preserving properties

[edgy sage]

@near trench So does $L_1 \leq_m L_2$ mean that if we have an algorithm that decides $L_2$, then there is an algorithm that decides $L_1$?

[pine jetty]

IlIIllIIIlllIIIIllll

[near trench]

yes

and if L1 is undecidable, all reductions also are.

[edgy sage]

So if $L_1$ is undecidable, it is still possible that $L_1 \leq_m L_2$ for some $L_2$

[pine jetty]

IlIIllIIIlllIIIIllll

[edgy sage]

But of course $L_2$ can't be decidable in that case

[pine jetty]

IlIIllIIIlllIIIIllll

[near trench]

Yeah

[edgy sage]

Then if $L$ is decidable, wouldn't we have $L \leq_m L'$ for every language $L'$

[pine jetty]

IlIIllIIIlllIIIIllll

[near trench]

That's the kind of fact that I was wondering existed, but how could a decidable language be reduced to an undecidable one? I was moreso thinking that maybe "every decidable language can be reduced to any other decidable language" kinda idea

[edgy sage]

If the language $L$ is decidable, then you don't need even need to use the algorithm for $L'$ to decide $L$

[pine jetty]

IlIIllIIIlllIIIIllll

[near trench]

I'm not sure I follow. Like in terms of languages and strings. the reduction is a function $f$ that maps one set of strings to another such that $W \in L $ iff $ f(W) \in L_1$, and the presence of a valid reduction is denoted $L \leq_m L_1$.

[pine jetty]

darkninja175

[near trench]

for a string W constructed from the relevant alphabet

I don't even really know how it works if you reduce a decidable language to an undecidable one, does that even work ? 🤔

[edgy sage]

The intuitive argument is clear. If a language is decidable, it trivially reduces to any other language since you already have an algorithm for the decidable one.

What is the formal definition for reduction?

Does it mean there exists a function $f : L \to L_1$ such that $w \in L \iff f(w) \in L_1$?

[pine jetty]

IlIIllIIIlllIIIIllll

[near trench]

I have this proof $L$ is decidable iff $L \leq_m L(000^(11+111))$, and I need that one idea that ties together the implication $L$ is decidable $\implies L \leq_m L(000^(11+111))$ y'know? because the reverse implication is easy. I'm getting caught up on thinking that I'm reducing any arbitrary language to something that seems so finite/specific.

The formal definition of a reduction is as I typed, yeah basically what you have there:

Let $A$ and $B$ be languages, then we say $A$ "reduces" to $B$ (written $A \leq_m B$) if there exists a computable function $f$ (an algorithm that always terminates) such that for every string $W$, $W \in A$ if and only if $f(W) \in B$

[pine jetty]

darkninja175

[edgy sage]

@near trench I think the difficulty is formalizing the function $f$. We already know that $f(w)$ should just ignore $L'$.

[pine jetty]

IlIIllIIIlllIIIIllll

[edgy sage]

How about define $f(w) \in L'$ if $w \in L$ and $f(w) \notin L'$ if $w \notin L$.

[pine jetty]

IlIIllIIIlllIIIIllll

[edgy sage]

Well pick $x \in L'$, $y \notin L'$. Set $f(w) = x$ if $w \in L$ and $f(w) = y$ if $w \notin L$.

[pine jetty]

IlIIllIIIlllIIIIllll

[edgy sage]

It seems that $f$ is computable since $L$ is decidable

[pine jetty]

IlIIllIIIlllIIIIllll

[near trench]

Hmm, I've never seen a construction like that but I guess it works?

[edgy sage]

It works as long as $L'$ is not $\emptyset$ and is not the entire set of strings.

[pine jetty]

IlIIllIIIlllIIIIllll

[near trench]

I was thinking that maybe another angle (trying for maximum formality) is to say. "well language L is decidable, and thus countable (not sure if that's a true implication), and L(000*(11+111)) is also decidable (thus countable), so a bijection could be formed between the two. This bijection is a valid reduction"

something like that idk

[edgy sage]

decidable does not imply countable

$\Sigma^*$ is decidable

[pine jetty]

IlIIllIIIlllIIIIllll

[stiff stone]

Can anyone help me with this? I tried using Rolle's theorem but I am unsure of the next steps.

[prime kraken]

maybe taylor's theorem would be more helpful

[stiff stone]

Ah ok, I will do that.

Sorry, I am not understanding how to apply it. Would you be able to give me another hint?

@prime kraken

I have an additional question. How is this differentiation matrix calculated? Is it a simple Jacobian matrix? Have tried looking online but can't seem to find it anywhere

[prime kraken]

it's made up from a finite difference scheme

precisely derived via the taylor theorem 😛

[stiff stone]

I am sorry, this is super confusing to me. Do you have any resource that might explain step by step?

[prime kraken]

if you look up those two topics, you should come across a nice explanation.

[stiff stone]

So is it a unique matrix for different function or is there some similarity between different matrices?

I always see the same example, -2 on the diagonal and 1 next to the -2.

[prime kraken]

https://www.dam.brown.edu/people/alcyew/handouts/numdiff.pdf

read this

[stiff stone]

Funnily enough I stumbled on this exact pdf an hour ago; however I am not sure how the matrix is created.

[prime kraken]

?

that's identical to the following

[stiff stone]

To be honest I'm not sure how to convert that into a matrix

or use that to make a matrix

[prime kraken]

$\begin{bmatrix} 1 & -2 & 1 & 0 & \dots & 0 \end{bmatrix} \begin{bmatrix} f_1 \ f_2 \ f_3 \ \vdots \ f_n \end{bmatrix}$

[pine jetty]

Edd

[prime kraken]

if we simply let f_1 = f(x), f_2 = f(x + \Delta x), and so on

[stiff stone]

Ok thats more understandable.

[prime kraken]

it's the same thing

a sum of scaled terms is a linear combination, and can therefore be written as a matrix-vector product

[stiff stone]

Ah ok

[prime kraken]

to get all of the desired finite differences, the nonzero terms [1 -2 1] have to be moved to the left and right. one simply puts these shifted versions of the vector as rows of a matrix

i'd recommend you brush up on your linear algebra and calculus

[stiff stone]

Alright thanks!

[fathom rain]

the coefficients are for example given here https://en.wikipedia.org/wiki/Finite_difference_coefficient#Central_finite_difference (you see that both your examples are second order, that is error is O(h^2))

[stiff stone]

Nice, thanks @fathom rain !

[naive moat]

holy mackerel

I've been reading about multigrid approaches to learn it

what I keep seeing is Brandt invented this, Brandt pioneered that, Brandt suggested an awesome algorithm that no one has been enough of a boss to implement yet, etc. EDIT: i see the question emoji, I don't have a question, just wanted to share that Brandt is apparently the go-to for multigrid approaches. Tho I have Briggs' Multigrid Tutorial which is suuuuper pedogagical

[fathom rain]

"no one has been enough of a boss to implement yet" I have worked on an industrial multigrid CFD code, theory and analysis of AMG, and also space-time MG. What do you mean?

[naive moat]

well, it's been implemented now, so I more mean recently

quote inspiring my paraphrase

from a paper by Mark Adams, where he implements a parallel version that can take advantage of the segmental refinement

[fathom rain]

No idea what that is. I worked on an edge-based RANS solver with agglomerated multigrid so it is trivially parallelizable and local

[naive moat]

there are a lot of very cool advances in multigrid methods and implementations in the past decade or so is what I'm discovering

[fathom rain]

(the code used to develop and design https://en.wikipedia.org/wiki/Saab_JAS_39_Gripen)

[naive moat]

Segmental refinement appears to be a 'bottom up' approach to multigrid where you can avoid storing the fine-grid solution/residual (and only store the coarsest-grid information) at the cost of some accuracy that comes from extra interpolation steps on the full solution (rather than interpolating the correction for prolongation, the normal method)

[fathom rain]

I personally just occasionally do analysis of MG nowadays. I do not enjoy "software engineering".

I think there is quite a lot of work currently on MG for non-local problems.

(and also space-time methods)

So, what do you want to use it for? 🙂

[naive moat]

Ah, I occasionally like to dive into the rabbithole of numerical methods and other applied math topics to keep myself refreshed of my options. I do electronic structure calculations and sometimes software development. A reading group I'm in is currently trying to implement Density Functional Theory, so I was thinking of learning-by-doing by using multigrid techniques for part of our implementation. Not sure where yet. Possibly the Kohn-Sham eigenvalue solving step, I've found a few papers on multigrid-based eigenvalue solvers.

My actual research is in modeling "strongly correlated systems" and I suppose I'm trying to see if multigrid techniques may be very useful in this sliver of my field

[naive moat]

I like seeing simple-but-effective algorithms. For example I recently implementated an algebraic variation of the "tetrahedron method" (so-called in condensed matter communities) that is simple to derive and code, yet still basically state-of-the-art in terms of efficiency. (It's a fancy integration routine for integrands in N-dimensions with integrable singularities)

Multigrid seems like another one of those simple-but-effective algorithms. I also don't enjoy "software engineering" so I like the simple-to-state methods since that can mean they're simple to implement (so long as you're not going for some state-of-the-art hybrid adaptive method or something)

[fathom rain]

Self-promotion, if you want a very simple (just like 100 lines of code) eigenvalue solver for structured matrices that will beat anything out there in efficiency (and is very accurate, pretty much machine precision), check out my "matrix-less methods" papers on https://www.2pi.se/publications/ (here is a TR with matlab code in it http://www.it.uu.se/research/publications/reports/2021-007/ some of the other papers have Julia code)

I am currently developing the same techniques for eigenvectors.

If you want me to try it on your "Kohn-Sham eigenvalue solving step" (maybe works?), then just send a function matrix(n) and I can take a look

[naive moat]

sometime sure when my reading group gets to implementing these things (in a week or two)

[prime kraken]

damn you have some good stuff there

[fathom rain]

thnx 🙂

[chilly elm]

quadratic spline, 4 points

[main flint]

Wait, why are there three interior points? Since there are only 3 splines, the continuity and differentiability will be checked at 2 points right?

That will make it 8 conditions, leaving 1 degree of freedom.

[chilly elm]

yes you're right

one think i dont get though

my book says continuity of spline itself doesn't give conditions

[main flint]

Sorry, I didn't understand. What conditions?

[chilly elm]

im talking in terms of number of conditions and number of unknowns

when we're counting for splines

[main flint]

So, if there are n points of data. So n-1 splines, continuity gets checked only at n-2 points.

If. I add this to the n constraints S_i = y_i, we will get a total of 2(n-1) constraints...

The only spline constructible with this much info is a linear spline.

So, even for quadratic splines with 3(n-2) variables, only continuity isn't enough.

[brave crypt]

Does anyone have experience with PINNs (Physics informed neural networks)?

[stiff stone]

How do I calculate the rate of convergence of the finite difference fixed point iteration method?

[wide spear]

Taylor series expansion probably

[floral shuttle]

what's a good (easy to implement) time integration for a stiff pde when doing fourier spectral collocation to solve it?

[wide spear]

Is it linear in time?

If so, any implicit method should be fairly straightforward and handle stiffness relatively well

[floral shuttle]

ah sorry

[floral shuttle]

i missed your message

it is linear in time but implicit is a bit of a pain

I need explicit to keep NlogN

implicit method would have me needing to solve a crazy system of eq'ns

[wide spear]

You want something explicit but robust to stiffness?

[floral shuttle]

yes

[wide spear]

Is there any sort of symplectic structure

[floral shuttle]

in particular i'm solving the diffusion brusselator system in 1D on a bounded domain

[wide spear]

Have you tried RK45

[floral shuttle]

yep

blows up after 3 iterations

since I'm time stepping in fourier space

[wide spear]

Lol rip

[floral shuttle]

lol exactly

[wide spear]

What about multistep methods

[floral shuttle]

haven't tried

[wide spear]

Try AB3 I guess

[floral shuttle]

you think it would converge better than rk45?

[wide spear]

No clue

If that doesn’t work you’ll need an implicit method

So you’ll want to look at ways of quickly inverting your matrix

Is it sparse

Does it have structure

How much does it change each time

[floral shuttle]

it does

good idea

bleh

i guess i'll need to do it that way

i was really trying to avoid implicit hahaha

[wide spear]

Anyways I’m sleeping now, I can try to help more in the morning

[floral shuttle]

np, appreciate it

[floral shuttle]

sooooo

rk4 is fine

I had to transpose my fft matrices 🤦‍♂️

stupid syntax error wasted 5 hours

[tribal pulsar]

does this method only apply to the Lagrange's Interpolating polynomial? does it work with the polynomial produced by Newton's difference method?

[wide spear]

Yes it should apply as well

I think

[wide spear]

To elaborate a bit more, n+1 points uniquely determine a degree n polynomial so Lagrange and Newton give the same polynomial

[tribal pulsar]

Thank you very much

[lament hazel]

am I crazy? Finding the weights of the Legendre polynomials by hand is nuts, right?

my professor alluded that he wants us to find the 5th degree Gaussian Quadrature "as exact as possible". I found the roots, which was not exactly easy but it was doable, but I'm starting to try to find the weights using Lagrange Interpolation and it's becoming a mess very quickly

most of the time you just use the table... right? Or have I just hit my computational ceiling lol

[prime kraken]

yeah these are the kinds of things you do once in small scenarios to convince yourself, then you toss it at a computer and never do it yourself again

[lament hazel]

okay glad I"m not crazy. It was not going well trying to solve it lol

[mint hemlock]

Yeah, intro numerical analysis courses having problems doing this stuff by hands is a little bonkers

I remember manually doing these quadrature methods as well as having to manually do a 9 point stencil in my intro numerical methods to ODEs/PDEs course

[forest musk]

hello. i'm currently struggling with the Digamma functions a little. I'd like to approximate them on the complex plane, but i'm not too sure how to do it relatively efficiently with arbitrary precision. I tried a bunch of formulas, including the Dirchlet integral one, which turned out to not be efficient enough, the Knuth phi(p/q) formula:

although, this isn't too efficient for obvious reasons - the last sum isn't convergent and since i usually feed 10^x as the denominator, it takes ages.

i know how to compute the values of phi for integer arguments - we could use the reflection formula to bring values from (-inf, inf) to (0, inf), then use the phi(x+1)=phi(x)+1/x & phi(1)=-gamma identities.

note that i'm not particularly interested in any polynomial expansions whose precision can't be adjusted in the runtime, since i want arbitrary precision.

right now i'm considering to use the phi(x) ~ ln(x) + 1/2x + O(...) identity, with a fixup for O(...) because the results are too far apart for small values of x.

what should do it is phi(x) = ln(x) + 1/2x - sum n=1 to inf of bernoulli(2n)/(2nx^(2n))

but i'd like to ask just in case if there's anything better or smarter

depending on how well this formula performs on smaller arguments, i believe i could shift the input range arbitrarily using the recursive identity, although this doesn't solve the problem of computing phi(z) for complex z.

[wide spear]

Mathematica implements digamma with Euler-Maclaurin summation, functional equations, and recursion

[forest musk]

I have a mathematica copy but i have no clue where to look for digamma implementation there. As for the Euler-Maclaurin formula, i already use it.

or rather plan to use it.

[wide spear]

Mathematica is closed source so you won't be able to find the actual code

I don't think

[forest musk]

isn't mathematica implemented in Wolfram language

there should be plaintext sources somewhere unless it's implemented in the kernel

as for SciPy, it seems like it uses the same approach i'm going to try out now.

[wide spear]

https://github.com/scipy/scipy/blob/701ffcc8a6f04509d115aac5e5681c538b5265a2/scipy/special/_digamma.pxd

Ah ok that's what I was going to suggest

[forest musk]

this seems to use the formula which uses Riemann zeta, which i haven't implemented either just yet.

oh well. i'll try to figure it out using Euler-MacLaurin and report back shortly.

[forest musk]

hold on... this formula works only for reals

[forest musk]

i decided to plug a few numbers into the formula I got from boost docs

unfortunately it seems very off

wikipedia seems to quote it too - maybe i'm missing something obvious?

[wide spear]

Is 10 terms enough

[forest musk]

not sure - wikipedia gives something around that as an example

see https://en.wikipedia.org/wiki/Digamma_function#Asymptotic_expansion

[wide spear]

Yes, those are the first few terms

Probably not enough though

[forest musk]

since these tend to be more accurate for larger z, maybe i could use the recursive property to bump it up when the input is small, say, between 0 and 20

[lament hazel]

hello good morning!
So my professor gave us this proof for a Gaussian Quadrature.
This is the problem:

This is the solution

My question is, since we have 4 unknowns, wouldn't simpson's 1/3 rule be inappropriate? I figured we would need to use the 3/8th rule, but according to his calculations, we got the correct result

When I did this problem, I did the 3/8ths rule in the same way and got the correct result.

So, since we have 4 unknowns, would Simpsons' 1/3 not be accurate?

Or is it because our basis goes form x^0 to x^3?

[turbid jay]

hey, can I use the interior point method for a constrained optimisation problem with a twice differentiable quasi-convex objective function?

[mint hemlock]

what do you mean that you used Simpsons’ 1/3? Simpsons’ is a special case of gaussian quadrature so it makes sense why’d you see a connectiont

[mint hemlock]

If you know the domain that the objective function is convex, then I don’t necessarily see the problem, but you also run into the same issues as non-convex functions with the interior point method, it’s just that now you are guaranteed to be convex on some region

with the region being of form (-\infty,a) iirc

where a is some fixed real number

[simple vigil]

in a first semester course trying to implement a 5point stencil for the 2d poisson problem

of course, here is the poisson problem (waiting on a plot to demonstrate that something is very much awry)

I believe my issues lie somewhere in my formulation... in particular, f and g. The gist of my code for solving the system (ie, Mv = F + G) is

h = 1/8
Ω = build_grid() # building (0,1) x (0,1)

u_mat = make_exact(u, Ω) # [0,1] x [0,1] of given exact solution
u_int = interior(u_mat)
x1_int = interior(Ω[1]) # build_grid is essentially matlab meshgrid
x2_int = interior(Ω[2])
    
M = h^2 * make_M()    
    
f = negLu(x1_int,x2_int) # - \Delta  u applied to the grid interior
g = h^2 * make_G(N,u_mat) # boundary values
  # make_G is essentially just a vector sum pulling the each edge [2:end-1] from u_mat

sum = f .+ g
S = reshape(sum', N^2, 1)
v = (reshape(M \ S, N, N))'

u(x1,x2) is given as sin(x1) * cos (x1) + (x2)^2. so, f(x) = 2sin(2x1) - 2

clearly this is not correct

code aside, I wanted to make sure that I was correct in choice for f and for g. Am fairly confident otherwise... got current code working with a 0 boundary condition example

[pine jetty]

mmmeeeker

[edgy sage]

You can apply the rule to any functions to get the weights. I think the standard way is to apply it to the Legendre polynomials $p_0, p_1, p_2, p_3$. You already know that the exact values are $(p_0, p_0), 0, 0, 0$ respectivley. @lament hazel

[pine jetty]

IlIIllIIIlllIIIIllll
Compile Error! Click the reaction for more information.
(You may edit your message to recompile.)

[simple vigil]

NVM this was an issue with my coefficients

[crimson phoenix]

isn't h close to something like 1e-5? in that case your error comes from that

[simple vigil]

@crimson phoenix h = 1/8 for that one

can talk more in a second... got my method working somewhat, but error is significantly higher than desired and not varying with 1/h^2 as it should

[crimson phoenix]

huh, that's a big problem

from your 5 point stencil, can you rebuild your discretization matrix correctly?

[simple vigil]

not sure what you mean, but here's what my discretization of the Laplacian looks like..

function discreteL(N::Int64)    
    db = I(N) * (-4)
    db = db + diagm(1 => ones(N-1), -1 => ones(N-1))

    M = kron(I(N), db) + diagm(N => ones((N-1)*(N)), -(N) => ones((N-1)*(N)))
    return -1 * M
end

right after I take this, I multiply by the 1/h^2

[crimson phoenix]

in principle something like $\begin{bmatrix}0 &1&0\ 1 & -4 & 1 \ 0 & 1 & 0\end{bmatrix}$

[pine jetty]

Anarchy

[simple vigil]

hm

I have 4's completely along the diagonal

[crimson phoenix]

talking about the stencil here

right

[simple vigil]

yeah

[crimson phoenix]

wait I wrote my internship report about a similar problem

[simple vigil]

(also, multiplied by -1 cuz I was implementing this for the poisson problem)

$\begin{cases} -\Delta u = f & x \in \Omega \ u = g & x \in \partial \Omega \end{cases}$

[crimson phoenix]

I mean, I'm okay with the discretization matrix

[wide spear]

This looks fine

You need to end cases

[pine jetty]

mmmeeeker

[simple vigil]

ah.. was puzzled for a second there

no wonder

[tribal pulsar]

Anyone could help me to verify the answer to this question? The solution that I get is in blue

It is under the topic of numerical differentiation

[simple vigil]

I think axes are flipped on the Error heatmap....

[crimson phoenix]

btw your discretization matrix should be divided by $h^2$

[pine jetty]

Anarchy

[simple vigil]

yea i do that right after I call it in the solver function

hold on

let me take another look at f and g

I haven't touched this in like a day or so

[crimson phoenix]

you have to do that before feeding the discretization matrix into the solver

[simple vigil]

yeah

h = 1 / (N+1)
    Ω = build_grid(0, 1, 0, 1, h)
    
    u_mat = make_exact(u, Ω)
    u_int = interior(u_mat)
    x1_int = interior(Ω[1])
    x2_int = interior(Ω[2])
    
    M = discreteL(N) / h^2
    
    f = negLu(x1_int,x2_int)
    g = make_G(N,u_mat) / (h^2)
    
    sum = f + g
    S = reshape(sum', N^2, 1)
    
    v = (reshape(M \ S, N, N))'
    return x1_int, x2_int, v, u_int, error_max_norm(u_int,v)

[crimson phoenix]

you're using matlab?

[simple vigil]

Julia

[crimson phoenix]

right

in u_mat I suppose you're computing the exact solution
How do you handle the boundary condition inside your discretization matrix?

[simple vigil]

so my understanding is that G is there to add the boundary conditions

[crimson phoenix]

I disagree

[simple vigil]

so what I'm doing is essentially taking the vectors along each side of the exact solution (ie, the boundary points) and adding those

so it looks like a ring of whatever boundary node is missing from the linear system in M

ie corners are the sum of the two boundary nodes missing from those, and edges are the values in the one boundary node missing from each of those respectively

[crimson phoenix]

I have a Jupyter notebook where I compute the solution of such a poisson problem, do you want it?

[simple vigil]

that'd be great

I'm convinced it's something I'm missing numerically

and that error just accumulates in the simulation

[crimson phoenix]

you should be able to plot the profile of your discretiztion matrix btw

you should have a diagonal coefficient that doesn't change

and "boxes" of coefficients

This image represents well what it should look like for your kind of problem

quickly looking at your coefficients, it looks like it's the right matrix, even though I don't have the full picture.

[fathom rain]

Easiest way to do 2D Laplace is just do 1D and use kron to construct the 2D

using LinearAlgebra

function laplace_1d(n)
    diagm(0=>-2*ones(Int64,n),-1=>ones(Int64,n-1),1=>ones(Int64,n-1))
end

function main()
    n1=10
    n2=20
    L1=laplace_1d(n1)
    L2=laplace_1d(n2)
    L2D=kron(L1,I(n2))+kron(I(n1),L2)
end

and if you want other BC just modify L1 and L2

[simple vigil]

took a nap.. back at it now

[simple vigil]

going to try this right now

[simple vigil]

I'm not completely sure what's meant by this... I think I'm just misunderstanding the Poisson problem from the get-go

@ Anarchy had mentioned that I was not correct in thinking that G from Mv = F + G is for the BC

I think I sent that somewhere way earlier... hm

### Create the G matrix
function make_G(N::Int64, u_mat::AbstractMatrix)
    g = zeros(N,N)
    g[1,:] += u_mat[[1],:][2:end-1]
    g[N,:] += u_mat[[end],:][2:end-1]
    g[:,1] += u_mat[:,[1]][2:end-1]
    g[:,N] += u_mat[:,[N]][2:end-1]
    return g
end

[fathom rain]

you modify the topleft and bottomright corners of L1 and L2 to get Neumann and top right and bottom left to get periodic

[crimson phoenix]

Listen to him rather than me, he's much better at this stuff than I am

[simple vigil]

ah got it..

what if I've got Dirichlet BC?

[fathom rain]

if you change -2-> -1 in a corner of L1 or L2 you get a Neumann BC there

Dirichlet you dont do anything (put the BC on the RHS if it is not homogenous BC)

[simple vigil]

(now that I think about it... I suppose I didn't specify Dirichlet BC earlier. mb) so, that is the G then?

[simple vigil]

.

or at least

This is what I have in mind, but 2d

where the [a,0,...,0,b]^T is G

[fathom rain]

IIRC it is just kron(RHS1,ones(n2))+kron(ones(n1),RHS2) to get the 2D RHS

[simple vigil]

@fathom rain do you mind if I just post/ send in totality

I am at wit's end

it's maybe like 100 lines and everything should be obvious in intent. the discretizations are the same

ie, mine - yours = 0

[simple vigil]

this yields some dimension issues, but I also don't quite get the reasoning... at this point, just trying to move on

[fathom rain]

sure

[brave crypt]

hello everybody
I'm trying to understund how instance encoding of a problem affect computational complexity
is there anyone i can ask some tips?

[wide spear]

http://www.matrixcalculus.org/matrixCalculus

Hello everyone

Just a cool little tool that I've found

[hollow pagoda]

Hi, I can't figure out how to do the dual of a summation, can anyone guide me?

The primal computes the min cost for a flow of 1 through the graph btw

[mystic orbit]

I'm looking to make a non-piecewise function with control points that would allow for something as diagramed

for reference, this is the original graph

[neat lion]

Is it feasible to compute the set of minimal forbidden minors for graphs of genus 1

[hollow pagoda]

My exam is in 2h welp I hope there isn't duals like this

[lofty bone]

I am reposting here, a question, I posted in #advanced-pdes :
Because generally evry one works in L², H1 spaces in case of numerical methods for PDEs, I am quite curious to know if there is Numerical methods for more exotic function spaces, or if one can apply it to those more exotic function spaces ?
By exotic I mean L^p-based Sobolev spaces, p not equal to 2, or even Triebel-Lizorkin or Besov spaces.

Just for general culture, I'm absolutely not related with the field of numerical PDEs

[crimson phoenix]

I personally never heard of anything other than L²-Sobolev spaces, or just L² to numerically solve PDEs
I have a masters degree in applied math for reference

Do SDEs count?

[wide spear]

Anatole I will answer your question when I get out of bed

[plain peak]

Hey all! Im currently trying to solve a problem that seemed relatively innocent, but turned out to be far more complex than I wished for...
Essentially Im trying to limit the rotation of a bone in a rig of an armature, based on 4 angles. Inside of the local space of a bone, the bone can be represented as a vector starting at the origin and the magnitude being the length of the bone. The default orientation of the bone (as defined by the rig) can then be understood as the "forward" axis of the coordinate system. Based on that forward axis you can make a unit vector (forward vector) and construct 2 other vectors, such that all of them are orthogonal to each other, who then together form the basis of the three dimensional vector space that this vector sits in. Now given 4 angles which rotate this forward vector around the origin "up","down","left" and "right" respectively you get 4 points that all sit on the surface of the unit sphere whos center also sits at the origin.

Now the first half of the problem is connecting these 4 points, such that it forms a "smooth" and closed boundary, where the forward vector always sits "inside" the boundary.

The second half of the problem is: Given an arbitrary point on the surface of the sphere, find the closest point on that boundary. Closest is here defined as "spherically" closest. So its not the euclidean distance but the "arc distance" between two points on the surface of a sphere.

My idea is using Quaternion Slerp to construct that boundary in a "Bezier-Curve like fashion" mimicking the De-Casteljau Algorithm and then somehow finding a function that lets you calculate the distance between a given point and the curve based on the interpolation value t, which I would then have to find the derivative of and solve for 0...

Any tips, ideas, or suggestions?

[wide spear]

Ok so

There are several different general classes of numerical techniques for pdes

1. finite difference methods: with these, you don't really consider what function space you're working in at all

2. finite element methods: with these, you compute a variational formulation and do indeed get a choice of function space, but we generally take the basis functions to be continuous everywhere (polynomials, piecewise linear) and the domain of the problem is compact so it doesn't make that much of a difference to just work in L^2 I think

3. spectral methods are just finite element methods with a particular choice of basis

I guess the main point is that we're working in a compact space and we want the solutions to be continuous and not just continuous almost everywhere so the various L^p spaces are not too different

[crimson phoenix]

If you're searching for a L^p solution, you're basically searching for a super regular solution, basically

[lofty bone]

what ?

No

Lp solution are absolutely not "more" regular

or at least I don't get what you mean

And thanks Ange, so if I am right, the numerical solvability of some PDE on Lp follows the same line, if the basis for solution is cnsistent in the Lp framework, Like Fourier series, polynomials on bounded sets, right ?

[wide spear]

Yeah

Like when working with pdes theoretically, the difference between the L^p spaces matters more because you actually want to consider things like singularities and weak solutions

But when you work with pdes numerically, you want no singularities and solutions that are differentiable in the classical sense

[lofty bone]

This explains why there is not that much explanations about the functional framework used in numerical PDE's expositions.

Thanks again

[cosmic thicket]

Is this a good place to ask about numerical methods?

[wide spear]

Sure

[cosmic thicket]

great thanks

Could you tell me what is a graded vs non-graded mesh of a numerical method and are graded meshes non-uniform?

[wide spear]

Do you take partitions to be uniformly spaced

[cosmic thicket]

I am not sure, this seems to describe a uniformly spaced mesh to me, what do you think?

can't see anything that explicitly tells me the step size

[wide spear]

If the partition is uniform then the graded mesh will not be uniform

[cosmic thicket]

can you tell me what "graded" means it is the first time I have come across the term? The expression in (3.5) in this image seems to describe something akin to a step size that differs for each step, which would be non-uniform...

Thanks for your help, if you do have an exact definition of a graded mesh please send it my way 🙂

[wide spear]

My understanding of a graded mesh is that it's a particular type of non-uniform mesh where you get more points around some regions of interest

I suppose this is a bit vague

But I don't think you'll find an exact definition

[cosmic thicket]

Yes that was what I was beginning to think. It is underestimated the extend to which math is an oral tradition

[fathom rain]

what book is this from?

[cosmic thicket]

by Charles Wing HO Green

[cosmic thicket]

Numerical Methods for Caputo–Hadamard Fractional
Differential Equations with Graded and Non-Uniform Meshes
[17:43]
by Charles Wing HO Green

[cosmic thicket]

Hi, Can anyone tell me for the A and B terms what is meant by t_k+1?

[mint hemlock]

I believe that it’s a side effect of the given weight coefficients a_j,k+1 and b_j,k+1 in order to construct A_j and B_j

as you’re looking at j from j=0 to k

so it’s just the k+1th timestep

[cosmic thicket]

Yes that is what I mean basically, how do you compute the k+1 time step if?

[mint hemlock]

Are you aware of implicit methods? (I’m assuming this is quadrature because I see trapezoidal rule)

[cosmic thicket]

I understand the method it is just the k+1 the timestep doesn't make sense to me

[mint hemlock]

oh, you just give it a timestep at k+1, I think

that is one potential way to find t_k+1

[cosmic thicket]

How does k+1 relate to t_N or T?

*relate

[mint hemlock]

t_k+1=t_1,…,t_N

because k = 0,1,…,N-1

[cosmic thicket]

I'm thinking..

This the first section before that and it has k=0,1,2,.....,N

[mint hemlock]

yes so you’re deriving y_k there
so when you find y_k+1, you need to allow k=0,1,…,N-1 as if k=N, t_N+1 doesn’t make sense in this context as the partition is [a,T] but log(t_N+1) > log(T) as log(t_N)=log(T)

so I think the book just assumed that this was clear and went with it and (through abuse of notation) redefined k

[cosmic thicket]

Thanks for helping, I don't get it yet 😦

[mint hemlock]

okay, so y_k=y(t_k)
so what does it mean by y_k+1?
well if k=N, then y_k+1=y(t_N+1) but that makes no sense, so it’s essentially just a change of variables doing
k := k-1

but just using the same variables

no, when k=N-1, t_k+1=t_N

which is T

it’s just redefining k so it makes sense to say y_k+1

[cosmic thicket]

Okay it is just a re-indexing type of situation, I must say it is written in a strange manner

[mint hemlock]

yes

[cosmic thicket]

Thank you for persevering with me @mint hemlock

[mint hemlock]

they just built up a definition for the predictor-corrector method without being careful with the definition of k

but revised that by the end by reindexing

yeah ofc

[cosmic thicket]

Mathematicians are can be too casual sometimes! 😠 😆

[cosmic thicket]

Can two sinusodials be manipulated into being parallel?

[wide spear]

What

[cosmic thicket]

never I am using fractions of circles instead

*mind

[stiff stone]

Can anyone help me attempt this question? I'm stuck on (a)

I'm not so sure if I have to do it by hand or using Matlab. For context we did the forward Euler method and the Crank-Nicholson method for odes.

[orchid sequoia]

you should ask that to your professor lol, imo use matlab: i wouldn't expect to solve this by hand from T=0 to T=2000

[nocturne valley]

how did it change from step 1 to step 2? this is taylor expansion series

[wide spear]

Solve (1) for f''(x0)

[nocturne valley]

im slightly unsure how to do that

[wide spear]

$f(x_0-2h)+f(x_0+2h)=2f(x_0)+4f''(x_0)h^2+\frac23\qty(f^{(4)}(\xi_1)+f^{(4)}(\xi_2))h^4$ so $4f''(x_0)h^2=f(x_0-2h)+f(x_0+2h)-2f(x_0)-\frac23\qty(f^{(4)}(\xi_1)+f^{(4)}(\xi_2))h^4$ so you divide through by $4h^2$ to get the second line

[pine jetty]

陆景和

[nocturne valley]

thank you so much

[stiff stone]

Why does this not work when I take f to be 2 or more functions for x and y but works for 2 or more functions of x and for 1 function of x and y?

\

[nocturne valley]

what does the epsilon mean here?

[frigid ravine]

It's a positive small number representing the absolute difference from the true function f and your polynomial P

[nocturne valley]

thank you

[wide spear]

Lorne Malvo it sounds like you have a coding question so ask in #computing-software

[brave crypt]

Maybe a bit weird question

Building subroutines for fluid dynamics solver that will run in cycles for my graduate thesis project

I had run code in c with malloc/free

And used almost the same skeleton from c in cpp with main differences being new/delete, and stdin/out for filewriting
I started learning cpp recently

Runtimes were 0.26 vs 18.9 seconds

Question, I thought of starting switching from C to CPP(which I barely know), but looking at the performance I doubt now. C was giving me nightmares with wrong memory allocations, messing up with pointers, but in terms of HPC running faster than Usain Bolt.

Do you recommend writing in pluses for the future anyway? as for more job research opportunities etc.

[fathom rain]

You are including compilations in those timings as far as I can see. Why are you building from scratch and not use something like FEnICS, OpenFoam, Deal.II? I don't think a workplace cares if you used c or c++

[brave crypt]

Idk , we’re not allowed to use third party software, maximum is IMSL/NR alike numerical libraries

[fathom rain]

ok, so what is the aim of the project?

you really should time the two without compilations, and as far as you get the correct result pick the fastest

[brave crypt]

do some navier stokes simulation, currently doing poisson solver

[brave crypt]

i still doubt that run without compile will take close to same time, cuz when i do manual g++/gcc compilation is almost instant for both

[fathom rain]

What kind of discretization to you use?

[brave crypt]

for poisson? regular central difference scheme

staggerred grid rie chow interpolation

i wanted to to collocated chorin projection, but will end up with some semi implicit most likely

[prime kraken]

yes, exactly the same way as in y = a + bx

[true hearth]

Hi sorry if this is a poor question but could someone perhaps explain this piece of text? I've read it so many times and still cannot understand what is going on

[mint hemlock]

@true hearth What's your question exactly about understanding

In essence, the Gaussian quadrature is the n+1th (I think) order integration method. So it'll work exactly up until a n+1th order polynomial.

One good way to think about it is you're interpolating the function f(x) using lagrange polynomials such that they agree with f(x) at the point $x_i$.
Then when you're integrating, think of it more like you're integrating
[\int_{-1}^1 f(x)dx \approx \int_{-1}^1 Q(x) dx= \sum_{i=1}^n \left(\int_{-1}^1 L_i(x)dx\right)f(x_i) = \sum_{i=1}^n c_i f(x_i)]

[pine jetty]

kirby

[lethal grail]

How do I know to write the Jacobi method in MATLAB?

I have a non-linear equation I have to approximate with an iterative method, and I want to see how the Jacobi method would be implemented in MATLAB, and how to know to do it

So that I would know how to use other algorithms and methods in MATLAB

[mint hemlock]

this is more of a #computing-software related question, but for approximating a non-linear equation, i’m curious how you’d use the jacobi method bc you need to decompose a matrix into a diagonal, upper, and lower matrix for the solution Ax=b. the algorithm implementation from there is fairly straightforward

while convergence not reached do
    for i := 1 step until n do
\sigma =0
        for j := 1 step until n do
            if j ≠ i then
\sigma =\sigma +a_{ij}x_{j}^{(k)}
            end
        end
x_{i}^{(k+1)}={{\frac {1}{a_{ii}}}\left({b_{i}-\sigma }\right)}
    end
k=k+1
end```

this was the algorithm from wikipedia

so you only need your base matrix A

[prime kraken]

you can avoid one of the loops by noting that you can add the diagonal to both sides of the equation and instead do an elementwise division or multiply by the inverse of that diagonal matrix (same thing). aside from that, you'd arrive at the linear equation either by doing a first order approximation and then multiplying by the transposed jacobian, or doing a second order approximation. you then do the jacobi method on that.

[cosmic thicket]

Can any one tell me why using particular graded (Non-uniform) meshes in a numerical method can lead to optimal convergence of the method?

[mint hemlock]

We can agree that if the mesh has infinitesimally small steps, convergence is going to occur, but this is impractical.
When referring to optimal convergence, we want something that is “close enough” in finite time. Essentially in regions where the derivative is large, you want the mesh to be in small steps, whereas if the derivatives are small (or smooth) larger steps suffice because you don’t need to worry as much about error

This is why non-uniform meshes are nice in that it gives us some level of control to the functions we’re working with

[pearl echo]

Btw do u know about mesh locking?

[mint hemlock]

was this towards me or oswald

bc I do admit that my answer did kinda neglect to mention that in the sense of convergence (although it does indeed converge iirc, just to a (much) lower solution) but I suppose that it’s nice to mention because the problem may not be as simple as “make the mesh small for better convergence”

but do you have a specific question regarding that or did you just want to point out this clarification?

[pearl echo]

this. There is this thing called mesh locking which basically states that errors increase cumulatively after each node. so it's possible to get errors with really small mesh.

[mint hemlock]

yes indeed, I should’ve mentioned that in my original answer, thanks for pointing that out :)

It’s mostly an order problem right? Like the use of lower order methods tends to give rise to such issues? I’ve only read about it in passing while I was doing some reading in numerical methods for PDEs

[pearl echo]

yes, whether the errors are positive or negative they always reduce the y-displacements in your solution.

that's like the crazy thing

[mint hemlock]

yeah, from what I can tell it’s a really interesting phenomenon

[nocturne valley]

Hi guys, i dont understand and know the comparisons of midpoint rule, gaussian , trapezium, simpson, monto carlo

i kknow why we use them but not too srue when to use them

and what are the advantages and disadvantages between them

[nocturne valley]

i dont understand these explainations

[fathom rain]

what part is unclear?

[nocturne valley]

i dont understand when to use the different quadrature methods

i understand when to use monte carlo

just not the rest

i know why to use them too

[fathom rain]

so the first case, MC is just very slow convergence, so proimarly used for high dimensional problem. this is not high dimensional. Trapezoidal rule evaluates the function, here 2/sqrt(x) at the endpoints and one end point is 0 and you can not evaluate there since it is infinity. hence you choose one of the two others, gauss will be more accurate (but more evaluations)

[nocturne valley]

oh wow, that makes so much more sense, thank you for your help

[fathom rain]

in second case midpoint, trapezoidal and gauss all give exact but midpoint is the least number of evaluations

third is high dimensional so MC

and the fourth as written

[nocturne valley]

awesome, thank you very much

[prime kraken]

what would be a good choice of preconditioner for a matrix of the form $\begin{bmatrix} 1 & -\gamma \ -\gamma & 1 \end{bmatrix}$, where $\gamma$ is complex?

[pine jetty]

19eddy4

[prime kraken]

i don't necessarily need to preserve symmetry, since the size is small

[candid timber]

is gamma a matrix?

[prime kraken]

just a complex scalar

it's a 2x2 mat

[candid timber]

oh

then can you take the inverse explicitly?

inverse is the ideal preconditioner

[prime kraken]

i can take it explicitly, but the results aren't nice

maybe i'm looking in the wrong place

lemme show an example of what i have, with some pretty pictures

so i'm trying to generate some scalar fields using the helmholtz equation and a 2 scatterer system. and somewhere down the line, that matrix shows up in an inverse problem

so, on the left, the scalar field in blue, and what it should look like in orange

in the middle image, what the spectrum of the scalar field in blue looks like

and on the far right, the condition number of that 2x2 matrix for each of the frequency bins

i'm computing the frequency bins of the scalar field exactly from the result of having inverted the matrix by hand

but no dice

it's the foldy-lax equations, if that means anything to you

or alternatively, a discrete version of the fredholm integral equation of the 2nd kind, or also lippmann-schwinger

[nocturne valley]

how would i do this because i dont understand the qs

[nocturne valley]

what is the derivation of the richardson extrapolation formula ?

[prime kraken]

just for completeness, i found a way to get it to work. no amount of regularization and preconditioning helped, but since the gamma entries in my matrix came from someone that should approximately represent an integral, using a different approximation scheme for that integral worked

just a lowly 2d trapezoid rule to compute new values for gamma

[burnt reef]

composite simpson's rule is an O(h^4) method, meaning the error is roughly proportional to h^4

https://en.wikipedia.org/wiki/Simpson's_rule#Composite_Simpson's_Rule

so $\frac{\text{Error when using } h_1}{\text{Error when using } h_2} \approx \frac{h_1^4}{h_2^4}$

[pine jetty]

jacob

[burnt reef]

comes from the error bound on the bottom

[nocturne valley]

thanks

[toxic wagon]

hello! i need some help in solving the 1d heat equation with a constant source being added to the whole area using crank-nicholson.
i already have some code (kind of) working for the base case where the heat applied is equal to 0 but im not exactly sure where to add it
for calculating the right hand side of the equation, our professor did a dot product with the previous timeframe and I'm not exactly sure how that would work

calculation code is as follows:
np and sp are just aliases of numpy and scipy respectively
k = delta t
h = delta x
g0 and g1 are the boundary conditions for u(0,t) and u(1,t) respectively

F = k / (2.0 * (h**2))
r = np.ones(m) * F
A = sp.sparse.spdiags([-r, 1.0 + 2.0 * r, -r], [-1, 0, 1], m, m).tocsr()
B = sp.sparse.spdiags([r, 1.0 - 2.0 * r, r], [-1, 0, 1],  m, m).tocsr()

for n in range(N - 1):
  b = B.dot(U[n, :]) + (k * f) # The (k * F) bit is where I'd assume this would go
  b[0] += F * (g0(t[n]) + g0(t[n+1])) # apply boundary conditions
  b[-1] += F * (g1(t[n]) + g1(t[n+1]))
  U[n+1, :] = sp.sparse.linalg.spsolve(A, b)

[lethal grail]

The thing about the Jacobi method was because I had some slides on an iterative method to approximate an equation, but I did not know how to write those steps in MATLAB.
I only asked about the Jacobi method to know how to "translate" or "interpret" something that is written on paper into MATLAB.
I eventually figured out writing the thing after a few more careful re-reads of the slides I had, but still had no idea on the convergence thing.
I think there was a lamba (that is within (0,1)) in the convergence criterion and two norms that were compared (the norm on the left side is less than or equal to that lamba multiplied by the norm on the right).
Anyways, I guess writing it... Is simpler than I thought it would be 5 days ago, since I thound the norm function in MATLAB (called norm, alongside the one for absolute values that I already knew of abs).
I will ask in #computing-software next time, if I ever need other help with being explained to why certain stuff is used and/or if an approximation method is still up to date in 2020-2022 or if it was in the 2010s.

Because I have a feeling that compression algorithms like singular value decomposition (SVD) had been superseded long ago (at least since the mid 1980s) when it comes to compressing images or even video frames, while still keeping all/most relevant information in it.

(I also tried a live demo of such a compression algorithm in real time... And could not understand how matrices or matrix mathematics was doing the things it was doing to the pixels. At all.)

[prime kraken]

depends what you mean by "better"

SVD is about as inefficient as most other data-driven methods in the sense that, although you get away with keeping only very few singular values, you also need the corresponding bases to reconstruct the approximate image

whereas bases generated from samples of a parametric model require fewer parameters

they may or may not yield the same reconstruction accuracy, or simply require more rank 1 matrices for the same acc

but overall need fewer params

depending on how you deal with parameters in ML, it could fall under either of these categories

and one would anyway use some sort of encoding on the result, too

[lethal grail]

In... Slightly simpler terms?

[prime kraken]

do you know DCT compression, for example?

[lethal grail]

For instance, let us compare video codecs.
For the same bitrate, AV1 (AOMedia Video 1) gives a better looking image than VP9 (which gives a better image than MPEG-4 AVC (H.264) does). And I will show an image of that.

[lethal grail]

Know, as in just hear of it or actually understand everything about it?

No, today is the first time I hear of discrete cosine transforms.

[prime kraken]

i can't really comment on much else, then

all of these codecs are far more sophisticated

[lethal grail]

Okay

So, I guess I will just never know about how computers understand how to use matrix stuff on individual pixels

[lethal grail]

rank, the thing that is determined with determinants that are different than from zero?

[prime kraken]

singular values or eigenvalues, yes

[prime kraken]

they don't, the procedure is coded by a person 😛

but as i said, a codec is rather sophisticated and does lots of things

just DCT or SVD compression is not nearly enough

[lethal grail]

Okay

Those two, too, are written in some sort of more complex software, right?

[prime kraken]

idk what you mean by that

[lethal grail]

I will try finding a publically available SVD thing in a programming language

And try explaining how SVD was used to transform the pixels the way it did

https://github.com/rameshputalapattu/jupyterexplore/blob/master/jupyter_interactive_environment_exploration.ipynb
https://medium.com/@rameshputalapattu/jupyter-python-image-compression-and-svd-an-interactive-exploration-703c953e44f6
https://gist.github.com/BLAKEMARG/af5f85627e5d6e557b17d0eb3a3f05bc
https://sites.utexas.edu/margolis/2019/05/15/image-compression-with-singular-value-decomposition/

[prime kraken]

most programming languages have a built in SVD routine

as for the actual implementation, you can just google them. i think the givens rotation approach is more or less straightforward

if you're not linear algebra savvy though, that will mean nothing

[lethal grail]

By the way, this is the live demo I tried (which did not bother explaining what was happening behind the scenes when I was sliding the cursor).
https://demonstrations.wolfram.com/ImageCompressionViaTheSingularValueDecomposition/

[prime kraken]

because it should be straightforward if you know what the SVD is

you'll have to look at the linalg

the application is straightforward

[lethal grail]

The thing is, when I had taken linear algebra and analytic geometry stuff I had taken in university, I was not going through a very good period of life.
On top of that, they were pretty sucky and I was not really in the mood to learn any of that.

But I do remember of the Givens rotation, by its name at least (not necessarily from those, but from approximation algorithms during numerical analysis, like numerical calculus and numerical methods).

[lethal grail]

If you say so...

[prime kraken]

what i would say is the crux here is that the SVD is equivalent to representing a matrix of rank R with a sum of R rank 1 matrices, each one associated to a singular value

you achieve compression by discarding the k matrices associated to the k smallest singular values, and you end up with a rank R - k matrix that is the "closest" to the original in some sense

the SVD is a decomposition of the form USV^T, where U and V are unitary/orthonormal, and S only has nonzero entries along its main diagonal. the number of nonzero entries in it is equal to the rank of the original matrix

and you set the smallest elements on the matrix S to 0

and the actual decomposition into USV^T can be done with givens rotations, as an example

if you anyway have to transmit the matrix or US'V^T, though, there is no gain

which is what i meant by having to "choose a nice basis"

like some fixed U and V that approximately diagonalize several images at the same time

or otherwise make the images sparse

(this is what the DCT does: one fixed choice of U and V^T makes many images sparse, and then U and V^T need not be transmitted)

[toxic wagon]

sorry, im bumping my question again

[foggy garden]

https://mathworld.wolfram.com/RSAEncryption.html
Question about the RSA math here, in an realistic situation, how does the receiver know the value of d?

[dense echo]

Hopefully he doesn't.

(For signatures, that is).

For encryption, the receiver of the message knows d because it's his private key and he generated it himself, by inverting e modulo phi(n).

[foggy garden]

Ok ty

that makes sense, I missed the fact that there were 2 pairs of keys

[foggy garden]

Another question, when the sender of the message encrypts the original message, does he use his public key or the recipients public key to encrypt the message?

[dense echo]

The recipient's public key if he's after confidentiality.

[zinc thicket]

I'm given the integral $\mathrm{I}=\int_{0}^{1} \log (1+\operatorname{sin} x) \mathrm{d} x$. Through the formula of Gaussian quadrature for 3 points, I can find an approximation to this integral. The question I'm asked is to show that
$$
2 \sup _{x \in[0,1]}|\log (1+\operatorname{sin} x)-P(x)|
$$
is an upper bound for the error of the Gaussian quadrature, knowing that $P(x)$ is the polynomial with degree $\leq 5$ that approximates $\log (1+\operatorname{sin} x)$ through least squares.
I don't even know where to start... Should I try to compute $P(x)$ and find the supremum? How would I show that it is an upper bound for the error?

[pine jetty]

ImHackingXD

[mint hemlock]

@zinc thicket What is the order of Gaussian quadrature using 3 points?

Using that, you will see the highest degree polynomial such that int P(x) = Gaussian Quadrature of P(x)

[zinc thicket]

I used the roots of Legendre polynomials, so 2*3-1=5

[tall solar]

Anybody here know if all "combinatorial optimization" problems correspond to matroids?

[sweet ore]

how to make a high level description of a TM to show that this is enumerable?
(it's akin to Goldbach's conjecture where every even integer greater than two is a sum of two primes)

[dense echo]

It's even computable. Given x you can just loop through all y's smaller than x and check whether y and x-y are both primes.

[pliant kelp]

Can anyone explain if this is perfectly secret or not?

I think its not perfectly secret

[dense echo]

That is a One-Time Pad. It is perfectly secure in theory, under the assumptions that the key is actually truly random (not just pesudo-random), the key is at least as long as the message, and the key is used only once. As soon as the same key is used for more than one message, everything falls apart.

https://en.wikipedia.org/wiki/One-time_pad

[pliant kelp]

yeah i figured that out after thinking about how similar it is to the OTP

Thanky ou !

[neon snow]

Calculating the total energy of a system at the initial state, and then using numerical integration until time T, then reversing all velocities, and iterating until time T (back to initial state) and now calculating the total energy of the system
Is this a good way to measure how well the numerical integrator conserves energy?

[mint hemlock]

Can you be more specific, like what numerical integration method are you performing?
If I'm understanding this correctly you have the following idea
E_0 -> integrate velocities over time - > E_f -> reverse velocities and integrate over the same time interval* -> E_1
So through the 2nd and 4th steps, you want E_0 and E_1 to be the same. The error is going to be O(order of the method), so you can calculate the maximal error analytically, or you can run it like this and see what E_1 - E_0 is.
*: This can also be formulated as saying "Going backwards in time from T to 0"

I am just not 100% sure of the question that you're asking here

[neon snow]

Hi! Thanks for your answer. I'm using velocity verlet (2nd order) and Yoshida 4th order.

And yes that was my idea. I want to compare these two methods quantitatively in my report, when simulating the solar system

My plan was otherwise to do som kind off root-mean square deviation to quantify energy conservation, since the energy oscillates around some value

[mint hemlock]

I guess the two important assumptions here are

1. Energy is conserved (which should be true, idk how you're implementing the solar system models)
2. The error from the integration methods are indeed related to the energy (as in, what are you integrating? Ig if you're approximating position from velocity, what energy are you finding or whatever)

Idk if that makes sense, my thing is just that okay, there will be integration error, which will present itself as some 3rd or 5th order term, but what conclusions you can draw from that, idk

[neon snow]

How can I quantify whether or not this is increasing over time?

Is there some kind of measurement that is good here?

[prime kraken]

you can consider plotting the rolling average

[neon snow]

How do I choose the parameter?

[prime kraken]

that's a very good question with no good answer

unless you know the underlying model

is this like a sinusoid with a time dependent vertical offset?

or some kind of random process + nonzero mean

[neon snow]

It's the energy of a symplectic integrator of our solar system

So the energy fluctuates around some approximation of the real energy

[prime kraken]

can you should the amplitude spectrum?

[neon snow]

What do you mean?

[prime kraken]

abs(fft(quantity you plotted))

and/or zoom into some part of the signal

cuz just looking at the plot you sent is not very useful on its own

[neon snow]

Yeah I know lmao

[neon snow]

ok I did this

[prime kraken]

plot that

[neon snow]

xD

[prime kraken]

zoom in near 0

[neon snow]

How far

[prime kraken]

as much as needed to see something useful

and i mean zoom along the x axis

[neon snow]

I'm finding stuff like this

I'm not sure how to interpret this

[prime kraken]

ok, this means your signal is largely harmonic

composed of sines and cosines

since you have well defined peaks

it means also that it should have 0 mean, and the mean should stay at 0 all throughout

then the rolling average should be good. as for how many samples for it, idk, try like 10% with 50% overlap between windows or something

or i think matlab has a "short time fourier transform"

which should effectively be the same

[zinc thicket]

In this derivation of a numerical differentiation formula, in the end the authors say that we can use an "alternative method" to find a common ksi. Does anyone know what that that method is?

[tall solar]

Does anybody here know about discrete optimization?

[mint hemlock]

what's up @tall solar

[tall solar]

Do you happen to know the difference between integer and combinatorial optimization?

I get the sense that there is a lot of overlap but that there is some difference

[prime kraken]

in some cases they are the same

but for example, integer opt can include the search for a vector whose entries are ints, but are otherwise unconstrained

while combinatorial problems usually involve a finite set where options can be drawn from without replacement

depending on the constraints, the problem can be translated back and forth in some cases

[tall solar]

Is linear programming a combinatorial problem because we know the objective has it's extrema at the corners (finite)?

I see linear programming in combinatorial opt textbooks and I'm trying to figure out why.

[prime kraken]

hmmm

i would say normally no

there is no reason to solve the problem that way in the first place

[tall solar]

Hmmm I'm gonna try to formulate my question better lol

[mint hemlock]

There are linear programming approaches to solving combinatorial problems yes

but it may be ill-posed to call linear programming a combinatorial problem

I think that it's mentioned bc linear programming is a crucial problem/idea in optimization (discrete or otherwise) bc of its formulation, where it can be framed as many problems, from optimal transport to optimizing polyhedra

[prime kraken]

that'd be my take too

[tall solar]

Ooooooooooh okay thank you that makes sense

[mint hemlock]

yeah, and if you have any other questions, feel free to ask :)

[random flame]

anyone has an example of a digitally applied partial differential equation for numeric computation?

maybe a physical simulation

or heat exchange

[mint hemlock]

wdym digitally applied

[random flame]

like in a simulation algorithm of some sort, maybe in a physical simulation that is rendered using OpenGL or anything of the like

[dense echo]

Atmospheric models for weather forecasting?

[prime kraken]

do you wanna see an example or do you just wanna know of some examples

[bright palm]

@onyx roost in here i think?

[onyx roost]

Hi smart people! I need some help understanding how to solve a generalized and non-linear eigenvalue problem using a cholesky decomposition to preserve symmetry. I think I have a solid grasp up till where the problem turns into $(A\lambda^2 + B\lambda + C) \times x = 0$ where $C = L^{-1} A (L^{-1})^T$ and $B = LL^T$. I don't quite know how to solve this.

[pine jetty]

Foureyed Jimmy

[onyx roost]

oop and the x is not a cross product haha sorry about that!!

[random flame]

anything i can get i'll be happy with

[random flame]

that's one possibility

[fathom rain]

https://nep-pack.github.io/NonlinearEigenproblems.jl/

[onyx roost]

Thank you so much!! This is perfect!!

[prime kraken]

does this help?

that's an example of the wave equation in a medium with 2 point-like inhomogeneities

[random flame]

ooo yes, exactly

do you have the algorithm/platform it was made in/for?

[prime kraken]

i made it in python

discrete approx to the green's function sol of the 2d helmholtz eq for several frequency bins, then ifft

[random flame]

how did you render it? i'm a C++ programmer, but recently getting into using matplotlib in python and others in R

[prime kraken]

i saved several png files and then a separate tool to turn that into a gif lol

[random flame]

ohh

nice

[prime kraken]

but there's the ion option for matplotlib

i think opencv might have better tools for it, matplotlib is kinda bad

[random flame]

i see

i wanna try and simulate either a heatmap using the heat differential equation, or something along those lines, maybe a vector interpolation with different speeds and forces considering fluid pressure or wave simulation

need to see the structure of writing a program around a partial differential equation so i can get started

[prime kraken]

you can look up something like "time domain finite differences"

[cosmic thicket]

Can anyone tell me where this quadrature formula comes from the paper doesn't specify?

[blissful fulcrum]

 Consider the string w = 0^k1^k. Since w ∈ L, it must be accepted by M. We think of processing w by M as a traversal of the corresponding directed graph.

Since the number of symbols in w is the same as the number of transitions made by M upon processing w, and the number of states in M is less that |w|, some state qi must be revisited while processing w.```
can some1 explain to me why ```the number of states in M is less that |w|``` this is true?

[dense echo]

This looks like it's in the middle of a longer argument.

Presumably k was defined on one of the preceding lines?

[rapid ridge]

Im trying to understand the BB84 protocal, and i dont get the last step where A creates a subset of n bits and sends it to B

and then they announce the bits of a and a' which corrospond

im unsure of what it means by this

like how do they know which bits they need to get rid of

or if they announce the corrosponding bits does this not compromise a and a'

[whole widget]

Hello! any book recos for the following? Thanks!
Stochastic Calculus
Stochastic Optimization
(I also asked in the book recos channel so if this post is considered spamming ill take it down)

[mint hemlock]

what sort of disciplines are you interested in? AI/ML, pure?
For stochastic calculus, Karatzas and Shreve seems good

[whole widget]

for AI/ML yeah, but id also like to read on the pure stuff later on. thanks good sir, ill start with Shreve

[mint hemlock]

Yeah, perhaps starting with https://web.stanford.edu/~jduchi/PCMIConvex/Duchi16.pdf for intro isn’t too bad

[prime kraken]

oh, this is a pretty nice, self-contained read

[rapid ridge]

was my question asked in the wrong channel?

[prime kraken]

well, it seems to be something of an error correction code

this is about as good as it'll get on this server. maybe try in a computer science or electrical engineering server

[brave crypt]

Hi everyone. I just joined discord after a couple of months. I am currently a 1st year PhD student working on Information Retrieval and Recommender Systems. i am interested in learning machine learning theory (maths side of it). If anyone is interested in the same, I'd be happy to connect 🙂

[rapid ridge]

ah ok, masybe a computer science one then but im not sure many of them will be able to awnser it either because it is on my maths course

but ill try anyway ty

[rapid ridge]

are you on about the more theoretical side of the coding side

i can code

but would call myslef a strong coder, only java,R,matlab

but id be intrested

[rapid ridge]

i have a lot more maths experiance tho, im the 3rd year of my maths degree, doing 2/3 applies 1/3 pure

[whole widget]

thanks so much!!!

[mint hemlock]

If you need specific books once you get through that, ask me, finding stuff that balances theory and practice is difficult sometimes, but typically I like to delve into the continuous case before I start looking at discrete applications, but you may be different

[mint hemlock]

Yeah, I found it to be good when I was first starting my work with gradient-free methods to just get notation and ideas down bc it's short

[brave crypt]

I am at a mix of both. But it's more applied work

[brave crypt]

Oh, that's cool

[brave crypt]

I am very interested in pure maths, but I don't have the necessary background to pursue it

[twin widget]

I recently stumbled into the wonderful world of Software Defined Radio, and by extension signal theory... it sparked quite a bit of fascination with the associated mathematics and I've been teaching myself more and more. And it's also very hands-on / applications-oriented.

[rapid ridge]

i think im combonation we might, i do the pure to further my applied and i do a fair mix

of pure subjects to do to that

[static peak]

is here the channel for numerical analysis

[mint hemlock]

yes

[earnest thicket]

hello, I am considering asking some computational linear algebra (grad level) questions here in the near future, is this the proper channel?

[mint hemlock]

yes

[prime kraken]

that seems about right

[mint hemlock]

Oh, Edd, I was 50% right about the timing, svd takes a long time, but most of it was bc I didn’t realize my data frame slicing was O(n^3)

[prime kraken]

[mint hemlock]

I forgot iloc[] was O(n)

[abstract zephyr]

I made a short blog to teach myself Gaussian Quadrature. Are there any cool concepts I should add?
https://observablehq.com/@laotzunami/numeric-integration-quadrature?ui=classic

[mint hemlock]

I think it's coming together pretty well, the only comment I have is to perhaps limit the nodes/weights of the legendre polynomials to n=50 because it'll still have the same effect, but it feels a little loud right now because it's so large

https://ieeexplore.ieee.org/document/9166469 @prime kraken Here's a paper I found that creates a dynamic time warping-like method using optimal transport (p-wasserstein). I'll look more into it, but it seems like a really interesting idea for some time series analysis

[prime kraken]

yo this is pretty cool

in case it helps, the one i had mentioned before is this one (kinda old) https://ieeexplore.ieee.org/stamp/stamp.jsp?arnumber=4567674

it's exclusively one to one and is only wasserstein-like, but maybe you find it interesting

[mint hemlock]

I will check it out in a little bit, my institution isn't showing up so I can't access this paper atm, I'll try on the vpn later

[prime kraken]

this one is really meant more for sparse time series, but still

[mint hemlock]

Yeah, that makes sense, but sparse methods are especially useful when the data is sufficient

so it’s an interesting idea nonetheless looking at the abstract

[small onyx]

please tell which channel is for numerical methods and optimisation

and how do I ask something in help channel it shows I don't have permission...

<@&268886789983436800>

[broken spoke]

this channel is probably the best place for that

but you should have permission if you use a channel in the "Available" section

[small onyx]

in conjugate gradient method isn't t the variable for g ? why are we changing x then ??

[prime kraken]

what t?

[small onyx]

g(x + vt)

x is fixed and then we find derivative of g and get g(x+vt)=g(x)-<v,b - Ax)^2 / <v,Av> by putting the t we get after finding minima of g and putting back in g

[prime kraken]

what you described is a procedure to find a new value of x

in your first image, all they did was define argmin

the argument of g for which g takes it's smallest value

and they called it x*

this is what you are looking for

and one way to do it is as you described, by taking your current value of x, finding some good direction v to move toward, and take a step of size t in that direction

[small onyx]

this is where I am confused "new value of x" wasn't x fixed and t the variable ?

[prime kraken]

no

or well, you can see it that way if you want, but if g is a function of x, then x + tv is effectively a new x

the goal in gradient-based methods is to find some optimal x* for which a target function g is minimal, i.e. argmin g(x)

and conjugate gradient does this by solving several smaller problems of the form argmin_t g(x + tv)

this is equivalent to forming a sequence x_k = x_{k-1} + t_k v which converges to x*

[fathom rain]

god I hate editorialmanager.....

[worn night]

gherkin method

[fleet sedge]

Hoho this is really interesting
https://github.com/zhangzheng2035/taot/blob/master/taot.py

[prime kraken]

nice

just make sure the license doesn't screw you over

[fleet sedge]

no license means all rights reserved ;_;

But honestly if I understood it well enough I think I can rewrite the thing

I'll be doing some time-series things soon but all I know about those is... it gets complicated real fast

[mint hemlock]

Yeah, time series data is really cool, but it does get complicated pretty quickly, but that's what time series analysis is for :) It's a pretty neat subject if you haven't taken a class/read anything on it

[brave crypt]

Hey everyone I am sophomore from Maths and Computing department. I happen to have Appled Computational Methods Theory and Lab in this semester. Can you guys suggest any good resources: online lec videos would be cool for this course??

[scarlet spruce]

I am constructing C^{k-2} interpolants and while this is trivial for odd degrees k, it's a bit problematic for even degrees (as in, I have an extra degree of freedom, and I am unsure what would be a good way to specify it). I found out that my cubic interpolation corresponds to using a binary primary pseudo-spline of order (2, 1), so I was wondering what the quadratic analogue would be.

Let k be the degree of the desired interpolant and n be a number of points (x1,y1), ..., (xn,yn), such that xi<x_{i+1}. Then for odd degrees I construct the C^{k-2} interpolant by stitching together pieces, such that the piece in [xi, x_{i+1}] is taken from the polynomial interpolating (x_{i-(k+1)/2+1}, y_{i-(k+1)/2+1}), ..., (x_i, y_i), (x_{i+1},y_{i+1}),...(x_{i+(k+1)/2},y_{i+(k+1)/2}) .

This poses no problem since there are an equal number of points on the left and right of the interval: (k+1)/2

for even degrees this is tricky since I could either pick the interpolating polynomial with 1 more point on the right, or one more point on the left. I would like to avoid such unsymmetric behaviour and have both of those accounted for. That is, for quadratic I want interpolation through the x_i and x_{i+1}points and somehow both x_{i-1} and x_{i+2} need to be taken into account for the third condition defining the quadratic polynomial.

My guess is that this would correspond to something like a binary primary pseudo-spline of order (1.5,1) but I haven't seen such a thing.

[pliant kelp]

Anyone have any tips for understanding how to use reduction to show that a certain construction of G is a PRG?

[pliant kelp]

If I have a PRG G that is pseudorandom and I have another G' that isnt psuedorrandom

If I were to XOR the output from those two generators

do I get a random generator?

garbage XOR with non garbage = ?

Isnt that just the one time pad

[dense echo]

By "isn't pseudorandom" do you mean it is truly random, or that it doesn't even look random?

[proud lion]

$\sum _{k=1}^{n-1}\left(\frac{1}{k}\right) is O\left(log\left(n\right)\right)$

I'm trying to prove this

[pine jetty]

The-Elite

[mint hemlock]

@proud lion recall that $\int \frac{1}{x}dx=\ln{x}$

[pine jetty]

kirby

[proud lion]

yep! I have that written down

I just not sure how to proceed with all this information because I think ln(x) >= summation right ?

[mint hemlock]

The sum is also bounded by the integral so

ofc the definite integral from 1 to x-1

but what is the definition of O(g(x))?

[proud lion]

f(n) ≤ cg(n) for some positive c and k, and for all n>=k

[proud lion]

@mint hemlock

[mint hemlock]

yup

sum 1/k <= C int_1^(n+1) 1/x dx

so sum 1/k <= C log(n+1)

You can fill in the details, but that’s essentially the way you go about it

[proud lion]

but that's not true right ?

how is that sum <= the integral

[mint hemlock]

Okay yes, the sum is <= 1+ int … for n>2

I forgot to add the constant argument

[proud lion]

dang i need it to work for k>=2

[mint hemlock]

that’s where you can select the constant

[proud lion]

$\sum _{k=1}^{n-1}\left(\frac{1}{k}\right) is O\left(log\left(n\right)\right)$

[pine jetty]

The-Elite

[proud lion]

the problem is that my summation upperbound is n-1

[mint hemlock]

It’s hard to be precise on my phone 😅

[proud lion]

so if I do n=2, log(2-1) = 0

[mint hemlock]

oh

[proud lion]

sum 1/k >= 0

thus ^

yeah

[mint hemlock]

I thought it was n+1, why do you need n>=2

[proud lion]

it's just the requirements for the summation i am working with

those are the bounds

[mint hemlock]

you just need it to be bounded above by Clog(x)

something like this

err, make that 2ln(x)

but the point stands

the first few n don’t really matter

because the definition only needs n>=L for some L (I changed it to L because you used k in your definition)

[proud lion]

ah i see i get it

[mint hemlock]

yeah, O(n),o(n), and \Theta(n) are asymptotic growths as n->\infty

so what happens in finite n doesn’t really matter too much

[proud lion]

and it wouldn't fail a proof either ?

if small n's didn't work

I guess if you choose your n value correctly, it's fine ?

[mint hemlock]

yes

[proud lion]

I see

okay thanks a lot

[mint hemlock]

the other two that I described have slightly different definitions, but the proof strategy works the same

[proud lion]

little o and theta ?

[mint hemlock]

yes

[proud lion]

sorry but if I was using k in summation(1/k), would I switch (1/x)dx to (1/k)dk ?

@mint hemlock

Sorry last @

[mint hemlock]

No worries

You can change it to k, just make sure your explanation is clear

[proud lion]

I see

thanks I appreciate all the help

[mint hemlock]

Yes ofc

[proud lion]

is it cool if i post what I have and maybe you can check it out ?

I think I finished

[mint hemlock]

Yes

That's fine with me

[proud lion]

Thanks

[mint hemlock]

This feels correct to me, the important thing to mention is the N such that this is true

where it holds for all n >= N

[proud lion]

true that

is there a good way to figure out what N should be ?

like is there another proof for that

[mint hemlock]

Essentially, the idea is choosing an n such that f(x) = g(x), then choosing the smallest integer above that

where f(x) is the original function you're working with and g(x) is its growth (what's inside of O(g(x)))

[proud lion]

I see okay

[naive yarrow]

math

[random flame]

this construction of an argument/solution is quite... bad, ik it's based on graph theory, but it's focused on an algorithmic argument. can i get any suggestions on how to structure this argument better?/more formally? it just feels like there's a lot wrong with it, and could be improved

[scarlet spruce]

An image would probably make the explanation redundant

As would an informal explanation