je joins un lien sur un programme SageMath que j'ai fais pour trouver la decomposition en valeur singuliere Algebrique (SVD en Anglais) de certaines matrices en dimension max ≤ 6 (peut fonctionner parfois en dimensions plus grande. SageMath a parfois des soucis pour trouver les eigen vecteurs en valeurs algebrique du fait (je pense) que l'algebre est trop compliqué) .
j'ai joins l'avertissement suivant dans le code:
beware that this little program is imperfect, (matrix C max dim <=6 )\
it is not always able to give algebraic values for U,V.\n\
even though SageMath knows how to calculate the algebraic values of the polynomials, \n \
but the search for the Eigen vectors for the 2 matrices U, V sometimes fails! \n \
I think the Algebra is sometime too much complicated \n \
in this case use Sagemath's u, s, v = M.SVD () function, which gives you numerical u,s,v \n \
example C=matrix(CDF,[(-2, 8, -7), (6, -3, -8), (-1, 8, -3)]); M.SVD()
Si quelqu’un a une idée pour l’ameliorer je suis preneur, je l'ai fais avec mes petites connaissances en codage SageMath, il doit y avoir plus simple.
sur Sage Cell
[edited le 2020/03/12 pour code Python 3]
je joins le code aussi ici:
- Code: Tout sélectionner
import itertools
def checkUV(C,u,s,v) :
# C=U*S*V^T ?
if C== u*s*(v.transpose()) :
return True
else:
return False
def build(combinaison0,eigenVectorList):
index=0
combinaison=list(combinaison0)
#show(combinaison)
#m=matrix(SR,vector(SR,[eigenVectorList[index][combinaison[0]]]))
m=(eigenVectorList[index][combinaison[0]])
index+=1
for i in range(1,len(combinaison)) :
#show(matrix(SR,vector(SR,[eigenVectorList[index][combinaison[i]]])))
m=m.augment(eigenVectorList[index][combinaison[i]])
index+=1
return m
def getSimplestUSVmatrices(C,CCTeigenVectorNormList,CTCeigenVectorNormList,s) :
u=CCTeigenVectorNormList[0][0]
for i in range (1,len(CCTeigenVectorNormList)):
u=u.augment(CCTeigenVectorNormList[i][0])
v=CTCeigenVectorNormList[0][0]
for i in range (1,len(CTCeigenVectorNormList)):
v=v.augment(CTCeigenVectorNormList[i][0])
show("u simplest:\t",u) ;show("Sigma simplest :\t",Sigma) ;show("v simplest :\t",v) ;
show("u*Sigma*(v.transpose()) simplest :\t",u*s*(v.transpose()))
if checkUV(C,u,s,v) :
usvList.append([u,s,v])
return usvList,True
##################### U
zeroCombinaison=[0]*(len(CCTeigenVectorNormList))
u=(build((zeroCombinaison),CCTeigenVectorNormList))
##################### V
zeroCombinaison=[0]*(len(CTCeigenVectorNormList))
v=(build((zeroCombinaison),CTCeigenVectorNormList))
#show("U S V^T :\t",u*Sigma*(v.transpose()))
usvList=[]
if checkUV(C,u,s,v) :
usvList.append([u,s,v])
return usvList,True
else:
return usvList,False
def swapCol(m,c_1,c_2) :
c1t=m.column(c_1)
m[:,c_1]=m.column(c_2)
m[:,c_2]=c1t
return m
def insert_row(M,k,row):
return matrix(M.rows()[:k]+[row]+M.rows()[k:])
def insert_column(M,k,column):
Mt=M.transpose()
m=matrix(Mt.rows()[:k]+[column]+Mt.rows()[k:])
return m.transpose()
def normalizeVectorList(eigenVectorList) :
eigenVectorNormList=[]
for l in eigenVectorList :
lambdaMatrixList=[]
for v in l :
v=vector(v)
#show(v/v.norm())
if v.norm() != 0 :
#show("v :\t",v,"v.norm() :\t",v.norm())
lambdaMatrixList.append(matrix(v/v.norm()).transpose())
else :
lambdaMatrixList.append(matrix(v).transpose())
eigenVectorNormList.append(lambdaMatrixList)
return eigenVectorNormList
def bonnetmatrixColumns(M):
"return a list of all possible columns arrangements"
colsArL=Permutations((M.ncols())).list()
#show(colsArL)
MbonnetedList=[]
for p in colsArL :
Margt=matrix(SR,[])
index=1
for c in p :
#show("c-1 :\t",c-1)
col=(M.columns()[c-1])
Margt=insert_column(Margt,index,col)
index+=1
#show(Margt)
MbonnetedList.append(Margt)
return MbonnetedList
#MbonnetedList=bonnetmatrixColumns(C)
#show(MbonnetedList)
def roundMatrix(M,rounding):
Mr=[[0 for c in (0..M.ncols()-1)] for r in (0..M.nrows()-1)]
for r in (0..M.nrows()-1):
for c in (0..M.ncols()-1):
Mr[r][c]=(round(M[r,c],rounding))
return matrix(Mr)
##########################################################################
def SolveNUorUNeqM(N,M,nu) :
# solve matrix N*U=M, with U Unknowns matrix
nR=N.nrows()
nC=N.ncols()
mR=M.nrows()
mC=M.ncols()
# generate unknown matrix variable list
#################################
uR=N.ncols()
uC=M.ncols()
uL=[]
for r in [0..uR-1] :
for c in [0..uC-1] :
uL.append("u_%d%d"%(r,c))
# workAound bug sagemath 8.8
# https://ask.sagemath.org/question/47470/var-list-with-only-one-variable-error/
# uV=var(uL) does not work if len(uL)=1
uV=[var(uL[0])] if len(uL) == 1 else var(uL)
#show("uV : ",uV)
#show(" U row : ",uR, " U col : ",uC)
Ul=[[0 for c in [0..uC-1]] for r in [0..uR-1]]
index=0
for r in [0..uR-1] :
for c in [0..uC-1] :
Ul[r][c]=uV[index]
#print "r : " ,r," c : " ,c ," index : ", index ," uV[index] : ", uV[index]
index+=1
U=(matrix(Ul))
#show("U : ",U)
#show("N : ",N)
#show("M : ",M)
#
if nu :
NU=N*U
else :
NU=U*N
eqT=[]
for r in [0..mR-1] :
for c in [0..mC-1] :
eqT.append(NU[r][c]==M[r][c])
# adding suplementary constraints if needed
#eqT.append(u_30==0)
########## end adding suplementary constraints
#for eq in eqT :
# show(eq)
S=solve(eqT,uV)
if len(S)<>0 :
show("Solutions : ",S)
if len(S[0])==U.nrows() * U.ncols() :
#for s in S[0]:
# show(s)
UnumL=[]
# fill the matrix Unum with its numerical values
UnumL=[[0 for c in [0..uC-1]] for r in [0..uR-1]]
index=0
for r in [0..uR-1] :
for c in [0..uC-1] :
UnumL[r][c]=S[0][index].rhs()
index+=1
Unum=matrix(UnumL)
# verify
#show("U : ",U," Unum : ",Unum)
#show("N : ",N)
#show("M : ",M)
#show("N*U :",N*U)
return U,Unum,S
else :
print "number of solutions do not match number of unknowns!"
return U,U,U
else :
print "No solutions !"
return U,U,U
"solution variables list (replace solution (sols) internal free variables r_i by arbitrary value 1 or 0) \
here this function gives vectors solution for the nullspace matrix (A-lambda_i*I)*X=0 \
in other words, this function gives a solution to the equation (A-lambda_i*I)*X=0 other than the trivial solution \
X=zero vector, this solution is the eigen vector corresponding to the eigen value lambda_i"
def solutionMatricesList(sols,X,lambd):
show("lambd : \t",lambd)
XcList=[]
varList=[]
varStr=""
for eqs in sols :
for eq in eqs :
#show(" eq : ",eq," eq.rhs().variables() : ",eq.rhs().variables())
for v in (eq.rhs().variables()) :
if v not in varList :
varList.append(v)
show(" varList : ",varList,"\t ",len(varList),"\tvariables")
if len(varList)==0 :
return XcList,varList
varListValue=[0] * len(varList)
Xc=copy(X)
for i in range(0,len(varList)) :
Xc=(Xc.subs({varList[i]:0 for i in range(0,len(varList))}))
XcList.append(Xc)
varListValue[0] = 1
for j in range(0,len(varList)) :
Xc=copy(X)
#show(" varListValue : ",varListValue)
for i in range(0,len(varList)) :
Xc=(Xc.subs({varList[i]:varListValue[i] for i in range(0,len(varList))}))
varListValue.insert(0,varListValue.pop())
XcList.append(Xc)
return XcList,varList
def getEig(A_lambdaI_list,lambdaList) :
rC=A_lambdaI_list[0].nrows()
zeroL=[0]*rC
ZeroVector=matrix(QQ,zeroL).transpose() ; # rotate vector horizontal to vector vertical
nu=True
#print "remember that A_lambdaI_list is the list of (A-lambda_i*I) matrices for all the lambda_i"
#print "now we will find the non-trivial (not zero vectors) solution for (A-lambda_i*I)*X=0"
index=1
lambdaIndex=0
eigenAllLambdaVectorList=[]
for e in A_lambdaI_list:
# solve matrix N*U=M
U,Unum,S=SolveNUorUNeqM(e,ZeroVector,nu)
#show(" N : ",N," U : ", U, " Unum : ",Unum ," = M : ",M)
#if nu :
# show(" N*Unum : ",N*Unum)
#else :
# show(" Unum*U : ",Unum*N)
XcList,varList=solutionMatricesList(S,Unum,lambdaList[lambdaIndex])
#show(" XcList : ", XcList)
show(" varList : ", varList)
eigenLambdaVectorList=[]
if varList <>0 :
for ee in XcList :
if ee <> ZeroVector :# remove the trivial one
eigenLambdaVectorList.append(ee)
show("A-(",(lambdaList[lambdaIndex]),")\t I")
show(" :",e," eigenV_" + str(index)+" : ",ee,LatexExpr(r"\lambda_"+str(lambdaIndex)))
index+=1
eigenAllLambdaVectorList.append(eigenLambdaVectorList)
lambdaIndex+=1
return(eigenAllLambdaVectorList)
#######################################################################################################
#######################################################################################################
########################## SVD Prcocess start here ####################################################
show(LatexExpr("C_{r.c}= U_{r.r} \,S_{r.c}\, V_{c.c}^T"))
show( "very good video on SVD :")
print "https://www.youtube.com/watch?v=EokL7E6o1AE"
show(LatexExpr(r"\text{as lambda is reserved Python word we use Lambda :} \, \Lambda \,\text{in place of} \, \lambda "))
Lambda=var('Lambda')
print "https://ocw.mit.edu/courses/mathematics/18-06sc-linear-algebra-fall-2011/" +\
"positive-definite-matrices-and-applications/singular-value-decomposition/MIT18_06SCF11_Ses3.5sum.pdf"
# ok for this matrix
#C=matrix(SR,[[4,8],[3,6]])
print "https://ocw.mit.edu/courses/mathematics/18-06sc-linear-algebra-fall-2011/"+\
"positive-definite-matrices-and-applications/singular-value-decomposition/#?w=535"
# ok for this matrix
#C=matrix(SR,[[5,5],[-1,7]])
# ok for this matrix
print "https://en.wikipedia.org/w/index.php?title=Singular_value_decomposition&action=edit§ion=7"
#C=matrix(SR,[
# [ 1 , 0 , 0 , 0 , 2 ],\
# [ 0 , 0 , 3 , 0 , 0 ],\
# [ 0 , 0 , 0 , 0 , 0 ],\
# [ 0 , 2 , 0 , 0 , 0 ] \
# ])
#
# OK for this matrix
#C=matrix(SR,[[3,2,2],[2,3,2]])
# ok for this one
print "https://www.youtube.com/watch?v=Ls2TgGFfZnU"
C=matrix(SR,[[1,1],[0,1],[-1,1]])
# ok for this matrix (SVD final matrix )
#C=matrix(SR,[[3,2,2],[2,3,-2]])
# ok for this matrix
print "http://web.mit.edu/be.400/www/SVD/Singular_Value_Decomposition.htm"
#C=matrix(SR,[[2,4],[1,3],[0,0],[0,0]])
#
## ok for this matrix
#C=matrix(SR,[ \
#[2, -1, 0, 0, 1],\
#[1, 0 , 0, 0, 0],\
#[0, 1 , 0, 0, 0],\
#[0, 0 , 1, 0, 0],\
#[0, 0 , 0, 1, 0] \
# ])
#ok for this one
#C=matrix(SR,[\
#(2, -1, 0, 1),\
#(1, 0, 0, 0), \
#(0, 1, 0, 0), \
#(0, 0, 1, 0) \
# ])
#C=matrix(SR,[(2, -1, 1), (1, 0, 0), (0, 1, 0)])
# test with a random matrix
# ok for this one
#C=matrix(SR,[(-5, 5, -1), (0, 2, 7)])
show( "beware that this little program is imperfect, (matrix C max dim <=6 )\
it is not always able to give algebraic values for U,V.\n\
even though SageMath knows how to calculate the algebraic values of the polynomials, \n \
but the search for the Eigen vectors for the 2 matrices U, V sometimes fails! \n \
I think the Algebra is sometime too much complicated \n \
in this case use Sagemath's u, s, v = M.SVD () function, which gives you numerical u,s,v \n \
example C=matrix(CDF,[(-2, 8, -7), (6, -3, -8), (-1, 8, -3)]); M.SVD() ")
# KO for these one
#C=matrix(SR,[(-2, 8, -7), (6, -3, -8), (-1, 8, -3)])
#C=matrix(SR,[(1, -9, 3), (-7, -4, -9), (3, -2, 1), (7, 8, 6)])
spanRC=5
spanValues=9
mR=ZZ.random_element(2,spanRC)
mC=ZZ.random_element(2,spanRC)
Ml=[[0 for c in [1..mC]] for r in [1..mR]]
for r in range(0,mR):
for c in range(0,mC):
Ml[r][c]=ZZ.random_element(-spanValues,spanValues)
" uncomment this line if you want test with random matrix "
#C=matrix(SR,Ml)
show("C matrix :\t",C,"rank :\t",C.rank())
##################################################################################
show ("get (to compute U) Eigen vectors for \t",LatexExpr(r"C \, C^T"))
CCT=(C*(C.transpose()))
show("CCT matrix :\t",CCT,"rank :\t",CCT.rank())
rC=C.nrows() ; cC=C.ncols()
Icct=matrix.identity(SR,rC,rC)
eigenValuesCCTEqu=(CCT-Lambda*Icct).det()==0
show("determinant of (CCT-Lambda*Icct) :\t",(CCT-Lambda*Icct).det())
CCTsol=solve(eigenValuesCCTEqu,Lambda,multiplicities=True)
if CCTsol==[] :
show(" no solution for (CCT-Lambda*Icct).det()==0 ")
exit()
CCTlambdaList=[]
index=0
for l in CCTsol[0] :
mul=CCTsol[1][index]
while mul > 0:
#The singular values are always real numbers
CCTlambdaList.append(real_part(l.rhs()))
mul-=1
index+=1
CCTlambdaList.sort(reverse=True)
show("CCTlambdaList :\t",CCTlambdaList)
CCTmatricesList=[]
for l in CCTlambdaList :
CCTmatricesList.append(CCT-(l)*Icct)
show("CCTmatricesList :\n",CCTmatricesList)
CCTeigenVectorList=getEig(CCTmatricesList,CCTlambdaList)
#######################################################################################################
show ("get (to compute V) Eigen vectors for \t",LatexExpr(r"C^T \, C"))
CTC=((C.transpose())*C)
show("CTC matrix :\t",CTC,"rank :\t",CTC.rank())
rC=C.nrows() ; cC=C.ncols()
Ictc=matrix.identity(SR,cC,cC)
eigenValuesCTCEqu=(CTC-Lambda*Ictc).det()==0
show("determinant of (CTC-Lambda*Ictc) :\t",(CTC-Lambda*Ictc).det())
CTCsol=solve(eigenValuesCTCEqu,Lambda,multiplicities=True)
if CTCsol==[] :
show(" no solution for (CTC-Lambda*Icct).det()==0 ")
exit()
CTClambdaList=[]
index=0
for l in CTCsol[0] :
mul=CTCsol[1][index]
while mul > 0:
#The singular values are always real numbers
CTClambdaList.append(real_part(l.rhs()))
mul-=1
index+=1
CTClambdaList.sort(reverse=True)
show("CTClambdaList :\t",CTClambdaList)
CTCmatricesList=[]
for l in CTClambdaList :
CTCmatricesList.append(CTC-(l)*Ictc)
show("CTCmatricesList :\n",CTCmatricesList)
CTCeigenVectorList=getEig(CTCmatricesList,CTClambdaList)
CCTeigenVectorNormList=normalizeVectorList(CCTeigenVectorList)
CTCeigenVectorNormList=normalizeVectorList(CTCeigenVectorList)
##########################################
###################" diagonal matrix Sigma"
LambdaList=CTClambdaList
LambdaDiag = diagonal_matrix(SR,CTClambdaList)
show(" diagonal_matrix with values^2 :\t", LambdaDiag)
SigmaDim=matrix(SR,rC,cC)
Isigma=matrix(SR,rC,cC)
for r in range(SigmaDim.nrows()):
for c in range(SigmaDim.ncols()):
if r < cC :
if r==c :
SigmaDim[r,c]= sqrt(LambdaDiag[r,c])
Isigma[r,c]= 1
else:
SigmaDim[r,c]=0
Isigma[r,c]=0
Sigma=SigmaDim
zeroUL=[0]*cC
ZeroUVector=matrix(SR,zeroUL).transpose()
zeroVL=[0]*rC
ZeroVVector=matrix(SR,zeroVL).transpose()
##########################################
#show("CCT (U) Eigen Vectors : ")
#show(CCTeigenVectorList)
for l in CCTeigenVectorNormList:
if l==[] :
l.append(ZeroUVector)
for l in CTCeigenVectorNormList:
if l==[] :
l.append(ZeroVVector)
show("CCT (U) Eigen normalized Vectors : ")
show(CCTeigenVectorNormList)
#show("CTC (V) Eigen Vectors : ")
#show(CTCeigenVectorList)
show("CTC (V) Eigen normalized Vectors : ")
show(CTCeigenVectorNormList)
CTClambdaMaskList=[]
plus=0
for i in range(len(CTClambdaList)):
if CTClambdaList[i] not in CTClambdaList[0:i] :
plus=0
CTClambdaMaskList.append(plus)
else :
CTClambdaMaskList.append(plus+1)
plus+=1
V=CTCeigenVectorNormList[0][0]
for c in range(1,len(CTCeigenVectorNormList)):
V=V.augment(CTCeigenVectorNormList[c][CTClambdaMaskList[c]])
show("CTClambdaList :\t",CTClambdaList)
show("CTClambdaMaskList :\t",CTClambdaMaskList)
#partial build of U
Upartial=(C)*(V[:,0])/sqrt(CTClambdaList[0])
for c in range(1,V.ncols()) :
if CTClambdaList[c]!=0 :
Upartial=Upartial.augment((C)*(V[:,c])/sqrt(CTClambdaList[c]))
show( Upartial)
U=copy(Upartial)
if Upartial.ncols()< len(CCTeigenVectorNormList) :
for c in range(U.ncols(),len(CCTeigenVectorNormList)) :
U=U.augment(CCTeigenVectorNormList[c][0])
show("U:\t",U,"Sigma:\t",Sigma,"V:\t",V)
show("U S V^T :\t",U*Sigma*(V.transpose()))
show("U S V^T round :\t", roundMatrix(U*Sigma*(V.transpose()),7)," C :\t",C)
[edited, Oops il y a avait une erreur sur la 5em ligne en partant de la fin :
j'avais ecrit: for c in range(V.ncols(),len(CCTeigenVectorNormList)) :
alors qu'il fallait ecrire bien sur (pour completer U avec le null space) :
for c in range(U.ncols(),len(CCTeigenVectorNormList)) :
je vais modifier le lien sage cell
c'est fait le lien sageCell est modifié.
PS, cet erreur n'intervenait pas dans le ca de la matrice utilisée par JJtheTutor car on n'avait pas a completer U , mais plantait le programme quand on avait a completer U !
.j'avais supprimer real_part sur les lignes 404,405:
- Code: Tout sélectionner
#The singular values are always real numbers
CCTlambdaList.append(real_part(l.rhs()))
et 445,446
- Code: Tout sélectionner
#The singular values are always real numbers
CTClambdaList.append(real_part(l.rhs()))
on peut les remettre, il est inutile de se trimballer l'Algebre imaginaire nulle dans les calculs.
une autre erreur
j'avais inverser les tailles des vecteurs, (sert quand pas de solution Eigen vecteur, excepté le vecteur trivial zero.!il faut mettre ligne 487 a 491
- Code: Tout sélectionner
zeroUL=[0]*rC
ZeroUVector=matrix(SR,zeroUL).transpose()
zeroVL=[0]*cC
ZeroVVector=matrix(SR,zeroVL).transpose()
a mon avis il doit y avoir encore d'autres erreurs, si vous en trouvez , merci de me les signaler ici.
Au debut je ne comprenais pas comment faire, je calculais les vecteurs Eigen pour U et pour V
et je tentai de construire naivement U et V avec ces Eigen vecteurs, mais cela ne fonctionnait pas pour toutes les matrices, parfois les signes des Eigen vecteurs n’etaient pas bon, et quand les valeurs Eigen sont multiple cela ne fonctionnait pas,et c’est JjtheTutor qui m’a fait comprendre comment faire.
Description de l’Algorithme:
calculer les valeurs Eigen
construire la matrice diagonal Sigma en rangeant les valeurs par ordre decroissantes.
Calculer les Eigen vecteurs pour U
calculer les Eigen vecteurs pour V
construire V en rangeant simplement les vecteurs Eigen
dans le cas ou il existe plusieurs vecteurs dans chaque liste de CTCeigenVectorNormList alors prendre successivement les different vecteurs de l’espace nul.
Example : pour :
C=matrix(SR,[ \
[2, -1, 0, 0, 1],\
[1, 0 , 0, 0, 0],\
[0, 1 , 0, 0, 0],\
[0, 0 , 1, 0, 0],\
[0, 0 , 0, 1, 0] \
])
on obtient la liste des Eigen vecteurs suivante :
ca fait ch.. ces P.. de <br> le compilateur Latex de ce site ne supporte pas les retour chariots dans le code Latex !
avec 3 fois les memes vecteurs pour la valeurs Lambda=1, dans ce cas il suffit de prendre un des vecteurs differents du null space pour les 3 cas (pour que chaque vecteur soit orthogonaux entre eux).l’ordre des vecteurs du null space importe peu.
c’est ce qui est fait avec le mask :
[0 , 0, 1, 2, 0]
correspondant au valeurs Lambda :
[3/2*sqrt(5) + 7/2, 1, 1, 1, -3/2*sqrt(5) + 7/2]
Calculer partiellement U UpartialColumn0=(C)*(V[:,0])/sqrt(CTClambdaList[0]),
UpartialColumn1=(C)*(V[:,1])/sqrt(CTClambdaList[1])
jusqu’a la dimension c de V rappel :C_{r.c}= U_{r.r} \,S_{r.c}\, V_{c.c}^T
ensuite completer simplement avec les Eigen vecteurs de U trouvées precedemment.
bon je ne suis pas sur d'avoir été tres clair !
pour plus de precision regarder l'excellente video de JJtheTutor:
excellente video de JJtheTutor