-2

I am using gaussian elimination to invert a matrix. Basically trying to use it to find the descent direction because somehow I feel taking the inverse is very slow. Gaussian elimination does speed up finding the descent direction. However during the course of the elimination if I start with particular initial vector, I am seeing that the last diagonal element is going to zero (essentially that means whole last row goes to zero). If I chose other starting points I don't get the issue but on few particular initial vector I get the issue. Has anyone encountered it and has ways to go around it? Also is there a faster way to find inverse of the matrix so that I don't have to do all this.

def powell_hessian(x):
    g1_ph = np.array([2.0 + 120.0 * (x[0] - x[3]), 20.0, 0.0, -120 * (x[0] - x[3])])
    g2_ph = np.array([20.0, 200.0 + 12.0 * (x[1] - 2.0 * x[2]) ** 2.0, -24.0 * (x[1]-2.0 * x[2]) ** 2.0, 0.0])
    g3_ph = np.array([0.0, -24.0 * (x[1] - 2.0 * x[2]) ** 2.0, 10.0 + 48.0 *(x[1] - 2.0 * x[2]) ** 2.0, -10.0])
    g4_ph = np.array([-120.0 * (x[0] - x[3]) ** 2.0, 0.0, -10.0, 10.0 + 120.0 * (x[0] - x[3]) ** 2.0])
    return np.array([g1_ph, g2_ph, g3_ph, g4_ph])

def powell_grad(x):
    g1_p = 2.0 * (x[0] + 10 * x[1]) + 40.0 * (x[0] - x[3]) ** 3.0
    g2_p = 20.0 * (x[0] + 10 * x[1]) + 4.0 * (x[1] - 2 * x[2]) ** 3.0
    g3_p = 10.0 * (x[2] - x[3]) - 8.0 * (x[1] - 2.0 * x[2]) ** 3.0
    g4_p = -10.0 * (x[2] - x[3]) - 40.0 * (x[0] - x[3]) ** 3.0
    return np.array([g1_p, g2_p, g3_p, g4_p])

def gauss_elimination(a, b):
    n = len(b)
    for k in range(0, n - 1):
        for i in range(k + 1, n):
            if a[i, k] != 0.0:
                lam = a[i, k] / a[k, k]
                a[i] = a[i] - lam * a[k]
                b[i] = b[i] - lam * b[k]
                print(a, "\n")
    # Back Substitution Phase
    for k in range(n - 1, -1, -1):
        b[k] = (b[k] - dot(a[k, k + 1:n], b[k + 1:n])) / a[k, k]
        # print(dot(a[k, k + 1:n], b[k + 1:n]), "\n")
        print(b, "\n")
    return b

dk = gauss_elimination(powell_hessian(x1), powell_grad(x1))

x1 = np.array([1.0, -1.0, 0.0, 1.0]) # creates the issue by making last row zero
x1 = np.array([3.0, -1.0, 0.0, 1.0]) # doesn't create the issue but takes longer to find the minimum
Mxxx_PFx
  • 1

2 Answers2

0

In your gaussian_elimination function the test

    if a[i, k] != 0.0:

should be

   if a[k, k] != 0.0:

However there is a more fundamental problem, which is that Gaussian elimination does not work for all invertible matrices. Consider the matrix

M = ( 0 1 )
    ( 1 0 )

This is invertible -- it is its own inverse -- but to use Gaussian elimination you need to swap columns.

Unless you have a particular interest in matrix inversion, you'd be better to use a library for inversion rather than write your own. Alas, I am to ignorant of python to recommend what to use.

dmuir
  • 4,211
  • 2
  • 14
  • 12
0

If you go through the gaussian elimination process for parallel vectors - or vectors which have parallel components - you get zeros along the diagonal (assuming ALL elements below the diagonal are parallel [equivalent]). If there is a zero along the diagonal then the algorithm likely will swap this row to the bottom, move the next diagonal element, and try to eliminate all remaining column elements for that row that was just pushed to the bottom in the process of reducing remaining rows. Unless there are other zero diagonal entries, then it will be successful in eliminating all of the remaining elements in the vector that was swapped to the bottom and you'll be left with a row of all zeros at the bottom.

At this point, if any given row is all zeros AND the answer (or augmented column element) for that row is zero, then you have a linearly dependent system. For each row such as this arbitrarily set the diagonal element to 1 (if desired - dependence implies you can make it whatever suits your needs).

If the row of the matrix is zero and the answer (or augmented column element) for that row is non-zero then this implies the system is not solvable.

*Rather than attempting to make hard rules of a the very general tool that gaussian elimination is, consider this useful facet of gaussian elimination method; every zero diagonal element remaining after reduction represents a dependent variable. You should convince yourself of this for matrices of various dimensions. You'll find it true for NXN matrices. For W>L append zero rows to form NXN matrix then count zero diagonals. For L>W instead cut bottom rows of zeros off then count diagonal zero entries.

wayne
  • 11
  • 2