search

uwhpsc-2016 / homework1

Homework #1
1 stars 1 forks source link

Exercise 3: solve_triangular #25

Open alyfarahat opened 8 years ago

alyfarahat commented 8 years ago

Are we expected to use solve_triangular() for the Jacobi method too? Is it okay to use element-wise inversion since S is diagonal? I am asking since the header of exercise3.py says that we only need to use the mentioned imported functions from numpy

cswiercz commented 8 years ago

You could use scipy.linalg.solve_triangular but it is unnecessary. Think of what operations you need to perform in order to solve the system

Dx = b

where D is diagonal. Since D is diagonal, but still dense, numpy.linalg.inv may still equate to a full, dense matrix inversion which is prone to the performance and stability issues mentioned.

As for the exercise3.py module "header" you are welcome to import any additional functionality you think you might need. I tried to already import the functions that will prove to be the most helpful.

cswiercz commented 8 years ago

Today (Tue, 12 April) there was a related discussion in the chat room between myself and @jhyearsley .

anders34 commented 8 years ago

Why would you have us import solve_triangular if we're not supposed to use it? Maybe I'm forgetting things from linear algebra, or maybe I'm just totally missing something, but why wouldn't you use it?

ckchmod commented 8 years ago

I think what the problem is trying to get at is that when you take the inverse of a square Diagonal matrix, it's simply the inverse of Diagonal entries: http://s-mat-pcs.oulu.fi/~mpa/matreng/eem1_5-3.htm

Notice that this is a unique case, and it wouldn't work when we're taking the inverse of (D + U). Which is why we would want to use linalg.solve_triangular

cswiercz commented 8 years ago

@anders34 Let me clarify: solve_triangular is unnecessary for the Jacobi method where the linear system you need to solve is diagonal. Try working out the math to see what you need to do to find x in

Dx = b.

It is, however, useful for Gauss-Seidel. The reason why you would want to use solve_triangular instead of numpy.linalg.solve is that the latter uses a "general" algorithm. As you will discover in the next homework (posted tomorrow) you can make some optimizations in the triangular case, which is what the former does.