Gram-Schmidt with arbitrary inner product?

edit

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,

load(eigen);
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.