1

Using numpy if I had a system of equations: 

3x + 2y = 5
1x + 4y = 10

I could solve them with numpy.linalg.solve: 

a = [[3,2],[1,4]]
b = [5,10]
solution = numpy.linalg.solve(a,b)

Now, what if I had two matrices, A and B each of shape (100,100) and I wanted to solve an equation of the form: (A'A - yBB')x = 0

I'm unsure as to how I would set this up using numpy.linalg.solve

atomsmasher
  • 725
  • 8
  • 20
  • Presumably `x` must be a vector with length 100, and `y` is scalar. Is that correct? Or is `x` also a scalar, as supposed in valentin's answer? – Warren Weckesser Jul 25 '16 at 17:34

2 Answers2

3

This looks like you are trying to solve the generalized eigenvalue problem, with y being the unknown generalized eigenvalue λ and x being the corresponding generalized eigenvector. If that is what you want, you can use scipy.linalg.eig. Here's an example.

To keep the output readable, I'll use arrays with shape (2, 2).

In [91]: from scipy.linalg import eig

In [92]: A
Out[92]: 
array([[2, 3],
       [1, 1]])

In [93]: B
Out[93]: 
array([[0, 1],
       [3, 0]])

These are the matrices in the equation.

In [94]: a = A.T.dot(A)

In [95]: b = B.dot(B.T)

Solve the generalized eigenvalue problem:

In [96]: lam, v = eig(a, b)

These are the generalized eigenvalues (your y):

In [97]: lam
Out[97]: array([ 6.09287487+0.j,  0.01823624+0.j])

The columns of v are the generalized eigenvectors (your x):

In [98]: v
Out[98]: 
array([[ 0.98803087, -0.81473616],
       [ 0.1542563 ,  0.57983187]])

Verify the solution. Note that the results are on the order of 1e-16, i.e. numerically close to 0.

In [99]: (a - lam[0]*b).dot(v[:,0])
Out[99]: array([  2.22044605e-16+0.j,  -8.88178420e-16+0.j])

In [100]: (a - lam[1]*b).dot(v[:,1])
Out[100]: array([ 0.+0.j,  0.+0.j])
Warren Weckesser
  • 110,654
  • 19
  • 194
  • 214
1

This specific equation has an analytical solution:

  1. either x is 0 and y can be anything

either

  1. x can be anything and then A'A - yBB' = 0 which can be linear and solved using numpy.linalg.solve

can be linear because the dimensions does not add well

valentin
  • 3,498
  • 15
  • 23
  • could you possibly help clarify something for me? the paper I'm reading states that the equation can be solved by the SVD of A and B. Would that then just be `np.linalg.svd(A)` and `np.linalg.svd(B)`, and I just take `s` ? – atomsmasher Jul 25 '16 at 13:31
  • maybe this link give you an insight http://stackoverflow.com/questions/32711920/what-is-benefit-to-use-svd-for-solving-ax-b – valentin Jul 25 '16 at 13:41
  • If `x` is assumed to be a scalar, then this is fine. However, the equation looks like the generalized eigenvalue problem, where `x` is a vector. – Warren Weckesser Jul 25 '16 at 18:00