Python Forum
Conjugate Gradient having issues with defining A (function to solve [A]{x} = {b} ) - Printable Version

+- Python Forum (https://python-forum.io)
+-- Forum: Python Coding (https://python-forum.io/forum-7.html)
+--- Forum: General Coding Help (https://python-forum.io/forum-8.html)
+--- Thread: Conjugate Gradient having issues with defining A (function to solve [A]{x} = {b} ) (/thread-25618.html)



Conjugate Gradient having issues with defining A (function to solve [A]{x} = {b} ) - DimosG - Apr-05-2020

Hi guys I have a problem with defining A (its in line 33) in my Full matrix solution and i was wondering if you could help me .

class Vector(object):
    def __init__(self, data):
        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(self.data, other.data)])

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

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

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

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


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:
                break
        beta = (rkk * rkk) / (rk * rk)
        pkk = rkk + beta * pk
        rk = rkk
        pk = pkk
        xk = xkk
        return xkk
    
#exit()

class FullMatrix(Matrix):
    def __init__(self, rows, cols, data):
        super().__init__(rows, cols)
        self.data = 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] += self.data[i+j*self.rows]*v.data[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*v.data[j]
        return Vector(data)
#exit()

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 = Af.cg(b)
        tf = time.time() - t0
        t_full.append(tf)
        print("full   : time: {0:6.2f}s (x[-1] = {1:12.10g})".format(tf, xf.data[-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 = As.cg(b)
        ts = time.time() - t0
        t_sparse.append(ts)
        print("sparse : time: {0:6.2f}s (x[-1] = {1:12.10g})".format(ts, xs.data[-1]))
        #exit()

        # 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, b.data)
        tn = time.time() - t0
        print("numpy  : time: {0:6.2f}s (x[-1] = {1:12.10g})".format(tn, xn[-1]))
        print()
        
    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_xlabel('size')
    ax.set_ylabel('times [s]')
    ax.grid()
    ax.legend()

    plt.show()



RE: Conjugate Gradient having issues with defining A (function to solve [A]{x} = {b} ) - micseydel - Apr-07-2020

Ideally you would be more specific about what you're looking for, and simplify your code so that it's only relevant to your question.


RE: Conjugate Gradient having issues with defining A (function to solve [A]{x} = {b} ) - 1968Edwards - Sep-18-2021

# Who solve sparse matrix with A sparce: A.x = f

# 2x + 3y + 0 = 6
# 0 + 0 + 2z = 4
# 4x + 0 + 1z = 2


#This is an example of sparse matrix above

# So I'm going to name rows, columns and data

import numpy as np
from scipy.sparse import bsr_matrix

row_bsr = np.array([0, 0, 1, 2, 2])
col_bsr = np.array([0, 1, 2, 0, 2])
data_bsr = np.array([2, 3, 2, 4, 1])

A = bsr_matrix((data_bsr, (row_bsr, col_bsr)), shape = (3,3) )
print(A)

# how to declare "b" in this case
# if I print A result OK
print(A)

# but I need x = A⁻1 f
# how do i proceed to extract the x

#=> This is my question because the scipy documentation doesn't show how to solve this.

# have many way:
# 1. bsr_matrix: Block Sparse Row matrix
# 2. coo_matrix: COOrdinate format matrix
# 4. csc_matrix: Compressed Sparce Column matrix
# 5. dia_matrix: Sparce matrix with DAgoal storage
# 6. dok_matrix: Dictionary Of Key based sparce matrix
# 7. lil_matrixÇ Row-based linked list sparse matrix.
# How to solve this matrix to extract the x.
# Common array resolution I can solve in Python but sparse
# arrays can't find the proper codes. The example above is
#pretty simple but my problem is 29x29 with a lot of zeros
#in A and only two values in f.