Gram-Schmidt with arbitrary inner product?

asked 2017-12-12 03:48:50 -0500 This post is a wiki. Anyone with karma >750 is welcome to improve it.

Hi,

I'd like to find an orthonormal basis of polynomials using an arbitrary inner product. I know how to do this in Maxima:

ip(f, g) := integrate(f * g, x, -1, 1); /* for example */
y : gramschmidt([1, x, x^2], ip);

But is there a nice way to do this in sage?

Thanks! -Andrew

edit retag close merge delete

Sort by » oldest newest most voted

In sage, there is a method gram_schnidt associated to a matrix.

After the matrix is declared, the method can be immediately called. The following defines a matrix A , computes the Gram-Schnidt matrices G.M as in the examples from

A.gram_schmidt?

and does a check for the relation $A=MG$, it is a possible starting point:

var( 'x' )
def ip(f, g):
return integrate( f * g, x, -1, 1 )

f_list = [1, x, x^2, x^3, x^4]
A = matrix( QQ, len(f_list),
[ ip( f, g ) for f in f_list for g in f_list ] )

print "A is\n%s\n" % A

G, M = A.gram_schmidt()

print "G is\n%s\n" % G
print "M is\n%s\n" % M

print "Is MG = A? %s" % bool( M*G == A )

Results:

A is
[  2   0 2/3   0 2/5]
[  0 2/3   0 2/5   0]
[2/3   0 2/5   0 2/7]
[  0 2/5   0 2/7   0]
[2/5   0 2/7   0 2/9]

G is
[         2          0        2/3          0        2/5]
[         0        2/3          0        2/5          0]
[ -424/5439          0  1376/9065          0   248/1813]
[         0    -12/595          0      4/119          0]
[ 64/112035          0 -128/22407          0    64/9603]

M is
[       1        0        0        0        0]
[       0        1        0        0        0]
[675/1813        0        1        0        0]
[       0   75/119        0        1        0]
[425/1813        0 950/1067        0        1]

Is MG = A? True

Later edit:

I try now hard to reproduce in sagemath the following result from maxima, obtained by copy+pasting into the maxima interpreter the code, thanks go to danielvolinski for this sample,

ip(f,g) := integrate(f*g, u, 0, 1);
y : gramschmidt([1, u, u^2, u^3, u^4], ip);

which delivers

2                           2
2 u - 1  6 u  - 6 u + 1  (2 u - 1) (10 u  - 10 u + 1)
(%o14) [1, -------, --------------, ----------------------------,
2           6                      20
4        3       2
70 u  - 140 u  + 90 u  - 20 u + 1
---------------------------------]
70

The situation is different, because there is no a priori knowledge of an orthonormal basis of the vector space generated by $1,u,u^2,u^3,u^4$ with the given scalar product $\langle f,g\rangle =\int_0^1f\;g$ .

To get a similar result i prefer to program from the scratch, there is no implemented facility...:

var( 'u' )
def ip(f, g):
return integrate( f * g, u, 0, 1 )

n = 5
B = vector( [u^k for k in range(n)] )
T = identity_matrix( QQ, n )    # in T we will have finally the transfromation matrix

C = vector( [u^k for k in range(n)] )    # in C we will have finally the Gram-Schmidt basis
for k in range(n-1):
S = identity_matrix( QQ, n )
for j in range(k+1, n):
S[j, k] = -ip( C[j], C[k] ) / ip( C[k], C[k] )

T = S*T
C = S*C

print T
print C

This gives the transformation matrix, and the new basis

[    1     0     0     0     0]
[ -1/2     1     0     0     0]
[  1/6    -1     1     0     0]
[-1/20   3/5  -3/2     1     0]
[ 1/70  -2/7   9/7    -2     1]
(1, u - 1/2, u^2 - u + 1/6, u^3 - 3/2*u^2 + 3/5*u - 1/20, u^4 - 2*u^3 + 9/7*u^2 - 2/7*u + 1/70)

in concordance with the maxima result.

more

Where is the resulting orthonormal basis of polynomials the user requested?

There is it in a simpler example:

sage: A = matrix( QQ, 2,2, [4,3,3,4] )
sage: G, M = A.gram_schmidt()
sage: A, G, M, M.inverse()
(
[4 3]  [     4      3]  [    1     0]  [     1      0]
[3 4], [-21/25  28/25], [24/25     1], [-24/25      1]
)
sage: -24/25*vector([4,3]) + vector([3,4])
(-21/25, 28/25)

Gram-Schmidt is applied on rows. The first row of A is copied in G. The lower triangular M gives the transformation in one direction. To obtain the second row in G, we use the linear combination from the M.inverse().

By the way, the user did not request any transformation, only a way to "do something similar", knowing that he knows how to do it in maxima. Doing it in maxima as he did, we also have to search for stuff, including transformation.

1

Hi Dan,

I think Andrew was looking for something similar to the following result in Maxima

Evaluate the following code in SageMath:

%maxima ip(f,g):=integrate(f*g,u,0,1)
%maxima y:gramschmidt([1,u,u^2,u^3,u^4],ip)

You will get five polynomials which are orthogonal (not orthonormal).

Is there something similar in SageMatt?

Thanks,

Daniel