I am trying to solve a matrix equation using a simple Minimum Residual Richardson method (http://depts.washington.edu/ph506/Boyd.pdf : page 304-306). I am using a finite difference matrix as preconditioner.

My initial matrix equation is: $$ Lu = F$$

where $L$ is not positive definite. To make it positive definite, I define $L_2 = L.L^T$. I use $L_2$ in place of $L$ hereafter. My finite difference preconditioner matrix is $H$. My initial guess $u_o = F$, the right hand side of the equation.

The algorithm I use is as follows:

$$ \begin{align} Lu &= F \\ Au &= g;\ where\ A = H^{-1}L,\ g = H^{-1}F\\ Residual:\ r_o &= g - Au_o \\ z_o &= H^{-1}r_o \\ \\ for\ k &= 1:max\_iter \\ \tau_n &= (r_n.A\ z_n)/(A\ z_n.A\ z_n) \\ u_{n+1} &= u_n + \tau_nz_n \\ r_{n+1} &= r_n - \tau_nA\ z_n \\ z_{n+1} &= H^{-1}r_{n+1} \end{align} $$

Using this algorithm, I get a solution which has correct shape but very high values. Is there some other method for making the matrix positive definite? Is my algorithm wrong? Or is my implementation incorrect?

This is what I get using my MRR method with finite difference preconditioning:enter image description here

I am implementing these on python. Here is a code snippet of what I am doing:

# Given the variable x, the preconditioner H, the matrix equation
# values L, F:

k = 30 # max iterations

L2 = np.dot(L, L.T) # to make the matrix positive definite

# New equation Au = g from Lu = f
# A = H^-1*L; g = H^-1*F
A = np.dot(np.linalg.inv(H),L2)
g = np.dot(np.linalg.inv(H), F)

# solution initial guess
uo = g
ui = g

# residual
ro = F2 - np.dot(A, uo)
ri = np.zeros((N))

res = np.zeros((k))

zo = np.dot(np.linalg.inv(H), ro)
zi = np.zeros((N))

for i in range(k):

    lamb = np.dot(A, zo) # temporary
    tau = np.dot(ro, lamb)/np.dot(lamb, lamb)

    ui = uo + tau*zo
    ri = ro - tau*np.dot(A, zo)

    zi = np.dot(np.linalg.inv(H), ri)

    res[i] = np.linalg.norm(ri)

    # update solutions
    uo = ui
    ro = ri
    zo = zi

For the entire code, you can check: https://github.com/prithvi-thakur/Jupyter_Notebooks/blob/master/MRR.ipynb


Your Answer

By clicking "Post Your Answer", you acknowledge that you have read our updated terms of service, privacy policy and cookie policy, and that your continued use of the website is subject to these policies.

Browse other questions tagged or ask your own question.