Python Scipy Error

Question

import scipy.sparse.linalg as scial
import scipy.sparse as scisp
import numpy



def buildB(A,x,col_size_A):
    d = numpy.zeros(col_size_A)
    for index in xrange(col_size_A):
        d[index] = 2*x[index]-1
    tmp = scisp.spdiags(d,0,col_size_A,col_size_A)
    return scisp.bmat([[A],[tmp]])


def buildQ(l,row_size_A):
    q = numpy.zeros(row_size_A)
    for index in xrange(row_size_A):
        q[index] = 2*l[index]
    return scisp.spdiags(q,0,row_size_A,row_size_A)

def buildh(A,x,b,col_size_A):
    p = A.dot(x)
    p = numpy.subtract(p, b)
    quad = numpy.zeros(col_size_A)
    for index in xrange(col_size_A):
        quad[index] = x[index]*x[index]-x[index]
    return numpy.concatenate((p, quad))


def ini():
    A = numpy.array([[1,1],[1,-1]])
    b = [1, 0]
    c = [1, 1]
    col_size_A = 2
    row_size_A = 2
    main(A,b,c,col_size_A,row_size_A)

def main(A,b,c, col_size_A, row_size_A):
    x = numpy.zeros(col_size_A)
    l = numpy.zeros(row_size_A*2)
    eps = 10e-6
    k = 0
    while True:
        B = buildB(A,x,col_size_A)
        Q = buildQ(l[row_size_A/2:row_size_A+1], col_size_A)
        Bt = B.transpose()
        h = buildh(A,x,b,col_size_A)
        g = numpy.add(c,Bt.dot(l))
        F = numpy.concatenate((g, h))
        print "Iteration " + str(k),
        tol = numpy.amax(F)
        print "- Tol "+ str(tol)
        if tol < eps:
            print "Done"
            break
        tF = -numpy.concatenate((c, h))
        FGrad2 = scisp.csc_matrix(scisp.bmat([[Q,Bt],[B, None]]))
        print FGrad2
        print FGrad2.todense()
        print " "
        print tF
        xdelta = scial.spsolve(FGrad2,tF)
        print xdelta
        x = x + xdelta[0:col_size_A]
        l = x[col_size_A:]
        k = k + 1


if __name__ == "__main__":
    ini()

The output is:

(2, 0)  1.0
(3, 0)  1.0
(4, 0)  -1.0
(2, 1)  1.0
(3, 1)  -1.0
(5, 1)  -1.0
(0, 2)  1.0
(1, 2)  1.0
(0, 3)  1.0
(1, 3)  -1.0
(0, 4)  -1.0
(1, 5)  -1.0
[[ 0.  0.  1.  1. -1.  0.]
[ 0.  0.  1. -1.  0. -1.]
[ 1.  1.  0.  0.  0.  0.]
[ 1. -1.  0.  0.  0.  0.]
[-1.  0.  0.  0.  0.  0.]
[ 0. -1.  0.  0.  0.  0.]]


lda must be >= MAX(N,1): lda=2 N=3BLAS error: Parameter number 7 passed to cblas_dtrsv had an invalid value
[-1. -1.  1. -0. -0. -0.]

So FGrad2 seems to be a valid csc matrix and tF a valid numpy.array.

What is wrong with this code? I don't even know why the error is before the print of tF even so the error is behind at spsolve

Edit

Ok i fixed that, it is because the first guess for the parameters was wrong leading to a singular matrix, but suppling a valid guess for l, leads to wrong calculation of spsolve

as mentioned i labeled all output as you can see spsolve returns the wrong calculation. $FGrad2 * xdelta != tF$

#!/usr/bin/env python
# -*- coding: utf-8 -*-


import scipy.sparse.linalg as scial
import scipy.sparse as scisp
import numpy



def buildB(A,x,col_size_A):
    d = numpy.zeros(col_size_A)
    for index in xrange(col_size_A):
        d[index] = 2*x[index]-1
    tmp = scisp.spdiags(d,0,col_size_A,col_size_A)
    return scisp.bmat([[A],[tmp]])


def buildQ(l,row_size_A):
    q = numpy.zeros(row_size_A)
    for index in xrange(row_size_A):
        q[index] = 2*l[index]
    return scisp.spdiags(q,0,row_size_A,row_size_A)

def buildh(A,x,b,col_size_A):
    p = A.dot(x)
    p = numpy.subtract(p, b)
    quad = numpy.zeros(col_size_A)
    for index in xrange(col_size_A):
        quad[index] = x[index]*x[index]-x[index]
    return numpy.concatenate((p, quad))


def ini():
    A = numpy.array([[1,1],[1,0]])
    b = [1, 0]
    c = [1, 1]
    col_size_A = 2
    row_size_A = 2
    main(A,b,c,col_size_A,row_size_A)

def main(A,b,c, col_size_A, row_size_A):
    x = numpy.zeros(col_size_A)
    x[0] = 0
    x[1] = 1
    l = numpy.ones(row_size_A*2)
    eps = 10e-6
    k = 0
    while True:
        B = buildB(A,x,col_size_A)
        Q = buildQ(l[row_size_A:], col_size_A)
        Bt = B.transpose()
        h = buildh(A,x,b,col_size_A)
        g = numpy.add(c,Bt.dot(l))
        F = numpy.concatenate((g, h))
        print "Iteration " + str(k),
        tol = numpy.amax(numpy.absolute(F))
        print "- Tol "+ str(tol)
        if tol < eps:
            print "Done"
            print x
            break
        tF = -numpy.concatenate((c, h))
        FGrad2 = scisp.csc_matrix(scisp.bmat([[Q,Bt],[B, None]]))
        print "FGrad2"
        print FGrad2.todense()
        print "tF"
        print tF
        xdelta = scial.spsolve(FGrad2,tF)
        print "spsolution"
        print xdelta
        print ""
        x = x + xdelta[0:col_size_A]
        l = xdelta[col_size_A:]
        k = k + 1


if __name__ == "__main__":
    ini()

Output:

Iteration 0 - Tol 3.0
FGrad2
[[ 2.  0.  1.  1. -1.  0.]
 [ 0.  2.  1.  0.  0.  1.]
 [ 1.  1.  0.  0.  0.  0.]
 [ 1.  0.  0.  0.  0.  0.]
 [-1.  0.  0.  0.  0.  0.]
 [ 0.  1.  0.  0.  0.  0.]]
tF
[-1. -1. -0. -0. -0. -0.]
spsolution
[-1. -1. -0. -0. -0. -0.]

Answer 1

The error I get is:

  File "stack27538259.py", line 62, in main
    xdelta = scial.spsolve(FGrad2,tF)
  File "/usr/lib/python2.7/dist-packages/scipy/sparse/linalg/dsolve/linsolve.py", line 143, in spsolve
    b, flag, options=options)
RuntimeError: superlu failure (singular matrix?) at line 100 in file scipy/sparse/linalg/dsolve/SuperLU/SRC/dsnode_bmod.c

As xnx wrote, FGrad2 is singular.

np.linalg.det(FGrad2.todense())   # 0.0

(scipy version 0.14.0)

---

after the change I get:

/usr/lib/python2.7/dist-packages/scipy/sparse/linalg/dsolve/linsolve.py:145:     MatrixRankWarning: Matrix is exactly singular

and

spsolution
[ nan  nan  nan  nan  nan  nan]

and an infinite loop unless I add k counter and break.

Answer 2

I think this is failing for you because your matrix is singular. E.g. convert to dense and use the regular numpy.linalg.solve:

>>> xdelta = numpy.linalg.solve(FGrad2.todense(), tF)

...
    raise LinAlgError('Singular matrix')
numpy.linalg.linalg.LinAlgError: Singular matrix

Answer 3

Documentation for cblas_dtrsv may be found (here)

Accordingly,

• the routine solves a triangular system A*X = B (presumably)
• lda is the leading dimension of matrix B
• N is the order of the matrix A
• the error message says lda = 2 and N = 3 but lda must be >= MAX(N,1)

Perhaps this helps track down the problem.