Conjugate Gradient having issues with defining A (function to solve [A]{x} = {b} )

Hi guys I have a problem with defining A in my Matrix where i def cg(self, b, epsilon=1e-4): and i was wondering if you could help me .

import math
import time
import matplotlib.pyplot as plt
import numpy as np

class Vector(object):
   def __init__(self, data): = data
       self.rows = len(data)

   def __mul__(self, other):
       assert self.rows == other.rows
       return sum([vi * vj for vi, vj in zip(,])

   def __rmul__(self, a):
       return Vector([a * d for d in])
   def __add__(self, other):
       assert self.rows == other.rows
       return Vector([i + j for (i, j) in zip(,])

   def __sub__(self, other):
       assert self.rows == other.rows
       return Vector([i - j for (i, j) in zip(,])

   def norm(self):
       return math.sqrt(self * self)

   def __str__(self):
       return '{0}'.format(

class Matrix(object):
   def __init__(self, rows, cols):
       self.rows = rows
       self.cols = cols

   def cg(self, b, epsilon=1e-4):
       ... # [1]
       xk = Vector([0.] * b.rows)
       rk = b - A * xk
       pk = 1. * rk
       for k in range(b.rows):
           alpha = (rk * rk) / (pk * (A * pk))
           xkk = xk + alpha * pk
           rkk = rk - alpha * (A * pk)
           if rkk.norm() < epsilon:
       beta = (rkk * rkk) / (rk * rk)
       pkk = rkk + beta * pk
       rk = rkk
       pk = pkk
       xk = xkk
       return xkk

class FullMatrix(Matrix):
   def __init__(self, rows, cols, data):
       super().__init__(rows, cols) = data

   def __mul__(self, v):
       assert self.cols == v.rows
       data = [0.] * self.rows
       ... # [2]
       for i in range(self.rows):
           for j in range(self.cols):
               data[i] +=[i+j*self.rows]*[j]
       return Vector(data)

#exit ()
class SparseMatrix(Matrix):
   def __init__(self, rows, cols, triplets):
       super().__init__(rows, cols)
       self.triplets = triplets

   def __mul__(self, v):
       assert self.cols == v.rows
       data = [0.] * self.rows
       ... # [3]
       for (i,j,d) in self.triplets:
           data[i] += d*[j]
       return Vector(data)

if __name__ == "__main__":

   sizes = [2 ** i for i in (2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12)]
   t_full = []
   t_sparse = []
   for N in sizes:
       print('_' * 80)
       print("size   : {0}".format(N))
       # b - vector
       data = [0. for i in range(N)]
       data[-1] = 10
       b = Vector(data)

       # Af - full matrix
       data = [0.] * (N * N)
       for i in range(N - 1):
           data[i + N * i      ] = 10.
           data[i + 1 + N * i  ] = -0.1
           data[i + N * (i + 1)] = -0.1
       data[N - 1 + N * (N - 1)] = 10.
       Af = FullMatrix(N, N, data)
       # solve
       t0 = time.time()
       #~ xf = cg(Af, b)
       xf =
       tf = time.time() - t0
       print("full   : time: {0:6.2f}s (x[-1] = {1:12.10g})".format(tf,[-1]))

       # As - sparse
       triplets = []
       for i in range(N - 1):
           triplets.append((i,     i,      10.))
           triplets.append((i + 1, i,     -0.1))
           triplets.append((i,     i + 1, -0.1))
       triplets.append((N - 1, N - 1,      10.))
       As = SparseMatrix(N, N, triplets)
       # solve
       t0 = time.time()
       #~ xs = cg(As, b)
       xs =
       ts = time.time() - t0
       print("sparse : time: {0:6.2f}s (x[-1] = {1:12.10g})".format(ts,[-1]))

       # An - numpy
       An = np.zeros((N, N))
       for i in range(N - 1):
           An[i, i] = 10
           An[i + 1, i] = -0.1
           An[i, i + 1] = -0.1
       # solve
       An[N - 1, N - 1] = 10.
       t0 = time.time()
       xn = np.linalg.solve(An,
       tn = time.time() - t0
       print("numpy  : time: {0:6.2f}s (x[-1] = {1:12.10g})".format(tn, xn[-1]))
   fig = plt.figure()
   ax = fig.add_subplot(111)
   ax.plot(sizes, t_full, 'b-o', label='full')
   ax.plot(sizes, t_sparse, 'r-o', label='sparse')
   ax.set_ylabel('times [s]')

Welcome, dimos.

What problem are you facing? A is not initialised anywhere, but is referenced in cg().

Hi ! Yes I know I want to create an A (square numpy matrix representing the linear operator.) But I just dont know how to define it in my code .

Why are you recreating vector and matrix objects rather than using a built in object?

I think you probably want a linear solver object that takes in a matrix (more generally a linear operator) and a vector and returns a vector with the solution. That’s a much better way to organize this sort of code.

Side note: If you are looking to create an efficient implementation of CG in Python… its been done by a ton of packages, like PETSc. I’d look at how those packages organize their code. If you just want to understand CG, then you are burning a lot of effort recreating matrices and vectors. Numpy has sparse matrices.

Hi Bob . I understand thats not the smartest way to go but the thing is our assignments is to solve this linear with 3 ways Full Matrix , Sparce, and numpy and then we want to compare how much time did it take with each one . And it’s a bit too much for me cause I started taking classes with Python 2 weeks ago so I’m really a beginner.

Ok, unless you are specifically told that CG must be a method of the matrix, then I would make a separate CG object. That will make your life much easier.

Unfortunately we were told to define the conjugade gradient into the Matrix … It’s okay honestly all the students are having an issue with this project . Thanks a lot for trying to help.

Huh, that’s a bizarre way to organize the code and certainly not what Numpy does, but I’m only a computational scientist, what do I know : )

I am looking at a few docs to see if there’s an obvious way to do it. But it is really a nonstandard way to organize code, for many reasons.

What is the exact requirement?

Lets see if I can translate it correctly in English . Write a program in Python programming language in which to define the classes Vector, Matrix (hyperclass), FullMatrix (subclass of Matrix, save with column layout), SparseMatrix (subclass of Matrix, save with list of lists) and write the above algorithm (linear solving with conjugate gradient) two cases of registers (FullMatrix and SparseMatrix). The solution of the system is required in an N-time diagram. And N goes from 2 up to 2^12.

See, to me that does not seem to require that the CG algorithm must be a method of Matrix. It just seems that you need a CG function that takes in your Matrix and Vector classes.

Well actually he gave us the code with missing pieces and in the Matrix there is
def cg(self, b, epsilon=1e-4):
and we need to write it ourselves.

Ah. That’s really strange.

Can you use self.__mul__(xk) in place of A * xk and so on?

1 Like

I think it works !!! Man thank you so much I think I destroyed my eyes staring this code for 2 days !

1 Like