how to fix "only length-1 arrays can be converted to Python scalars" when getting integral with an array argument

Question

I am using quad from scipy.integrate to get an integral in a limited range from an object. suppose the target object is in the blow:

∫expm(A*X).expm(B*X)dx

which both A and B are numpy matrix.

To solve this I have used blow code:

from scipy.integrate import quad
from scipy.linalg import expm
import numpy as np

def integrand(X, A, B):
    return np.dot(expm(A*X),expm(B*X))


A = np.array([[1, 2], [3, 4]])
B = np.array([[1, 2], [3, 4]])

I= quad(integrand, 0, 1, args=(A,B))

But for the result I get this error:

TypeError: only length-1 arrays can be converted to Python scalars

I know that The error "only length-1 arrays can be converted to Python scalars" is raised when the function expects a single value but you pass an array instead. but my problem is based on array. so how can I fix it.

Double check the docs regarding the function argument to quad. — hpaulj, Feb 16 '19 at 16:14
Don't you know where to find the documentation for the Python code that you use? How did you even find out how to call `quad` as much as you did? — hpaulj, Feb 16 '19 at 16:53
OK, let's be more explicit - `quad` integrates a scalar function. Your function returns an array. — hpaulj, Feb 16 '19 at 16:56
I'm new in using 'quad' and I've found it by searching and then reading related document from 'scipy' website. but they are just simple explanation with one example. if i said could you recommend me, i mean more detailed documents. @hpaulj — samie, Feb 16 '19 at 17:42
Do you want separate one-dimensional integrations for each element of `A` (and `B`) or do you want to multiply a 2-D vector `x` with the matrix and integrate in 2 dimensions? The latter would be inconsistent, since you cannot add the vector `A*x^2` and the matrix `B`. — ml4294, Feb 17 '19 at 14:21

score 1 · Accepted Answer · answered Feb 18 '19 at 13:31

As pointed in the comments, quad expects a scalar function. You can always pass the function to a scalar by adding the index as an output:

def integrand(X, A, B, ix=None):
    """ pass ix=None to return the matrix, ix = 0,1,2,3 to return an element"""
    output = np.dot(expm(A*X),expm(B*X))
    if ix is None:
        return output
    i, j = ix//2, ix%2
    return output[i,j]
I= np.array([quad(integrand, 0, 1, args=(A,B, i))[0]
for i in range(4)]).reshape(2,2)
I
>>array([[1031.61668602, 1502.47836021],
       [2253.71754031, 3285.33422634]])

Note that this is very inefficient since you are calculating the integral 4 times, as long as this doesn't bother you.

Alternatively, use trapz:

x_i = np.linspace(0,1,60)
np.trapz([integrand(x, A, B) for x in x_i], x=x_i, axis=0)
>>array([[1034.46472361, 1506.62915374],
   [2259.94373062, 3294.40845422]])

yes I have used your idea to solve the problem by indexing but it is too slow for large matrix and more complex integration. So I tried your recommendation in alternative way and the speed is improved. it's really useful and great — samie, Feb 18 '19 at 16:36

Nico Schlömer · Answer 2 · 2020-02-03T12:52:23.963

quadpy does vectorized computation. The fact that expm only works on square matrices (and not on lists of square matrices) requires a bit of juggling with the matrix shapes, though.

from quadpy import quad
import numpy as np
from scipy.linalg import expm


A = np.array([[1, 2], [3, 4]])
B = np.array([[1, 2], [3, 4]])


def integrand(X):
    expAX = np.array([expm(A * x) for x in X])
    expAX = np.moveaxis(expAX, 0, -1)
    #
    expBX = np.array([expm(B * x) for x in X])
    expBX = np.moveaxis(expBX, 0, -1)
    return np.einsum("ij...,jk...->ik...", expAX, expBX)


val, err = quad(integrand, 0, 1)
print(val)

[[1031.61668602 1502.47836021]
 [2253.71754031 3285.33422633]]

how to fix "only length-1 arrays can be converted to Python scalars" when getting integral with an array argument

2 Answers2