MATH920 L09: Nonlinear DNN model reduction

Read data

SetDirectory[NotebookDirectory[] <> "/models/1d/case"];
file=FileNames["*.out"][[1]];
stream = OpenRead[file,BinaryFormat->True];
m=BinaryRead[stream,"Integer32"]; (* Nr. of point masses *)
n=BinaryRead[stream,"Integer32"]; (* Nr. of time steps *)
uT=Table[0,n,m];
For[k=1,k≤n,k++,
For[j=1,j≤m,j++,
uT[[k,j]]=BinaryRead[stream,"Real64"];
]
]
Close[stream]; u=Transpose[uT];
{{m,n},Dimensions[u]}

Manipulate[ListPlot[u[[;;,kt]],Frame->True,PlotRange->{{0,m},{-2,2}},GridLines->Automatic],{kt,1,n,1}]

Extract principal modes

Center the data

μ = Mean[#]& /@ u;
uC=u;
For[i=1,i≤m,i++,uC[[i]] = u[[i]] - μ[[i]]];

Dimensions[uC]

Table[Chop[Mean[uC[[i,;;]]]],{i,1,m}]

Dimensions[uC]

Carry out SVD of centered data to obtain principal components of correlation

{Q,Σ,R}=SingularValueDecomposition[uC];

Check that singular values indicate dominant modes

ListPlot[Log10[Diagonal[Σ]],PlotRange->All,Frame->True,GridLines->Automatic,PlotLabel->"Potential energy spectrum"]

Manipulate[ListPlot[Q[[;;,l]],Joined->True,Frame->True,PlotRange->{{1,m},All},GridLines->Automatic],{l,1,Min[m,n],1}]

DNN model reduction

Linear layers, is the SVD recovered?

DOF reduction at each stage

Define the number of degree of freedom reduction per layer δDOF,  final number of DOFs, thus defining the number of layers

δDOF=2; mFinal=6; nLayers=(m-mFinal)/δDOF;
lm = NetChain[
Table[LinearLayer[m - k δDOF,"Input"-> m - (k-1) δDOF,"Biases"->None],{k,1,nLayers}] ]

Define SVD reduced basis

QmFinal = Q[[;;,1;;mFinal]]; QmFinalT=Transpose[QmFinal]; Dimensions[QmFinal]

Define data set

Dimensions[u]

Dimensions[QmFinal]

Dimensions[outputs]

Dimensions[uT[[1]]]

inputs = uT;
outputs =  uT . QmFinal;
data = Table[inputs[[k]]->outputs[[k]],{k,1,n}];
Dimensions[#]& /@ {inputs,outputs}

Carry out SGD to find DNN weights

LinearModel = NetTrain[lm,data]

Define a function to extract the layer weights

W[l_]:=Normal[NetExtract[NetExtract[LinearModel,l],"Weights"]];
{Dimensions[W[1]],Dimensions[W[nLayers]]}

NNT=W[nLayers];
For[k=nLayers-1,k≥1,k--,NNT=NNT . W[k]];
NN=Transpose[NNT];
Dimensions[NN]

Does the trained DNN carry out projection onto the SVD dominant modes? First, visualize the DNN modes

Manipulate[ListPlot[NN[[;;,l]],Joined->True,Frame->True,PlotRange->{{1,m},All},GridLines->Automatic],{l,1,mFinal,1}]

No, from the above visualization the eigenmodes are not obtained. Orthonormalize the DNN linear transformation to check if dominant modes are obtained

{QNNT,RNN}=QRDecomposition[NN]; QNN=Transpose[QNNT]; Dimensions[QNN]

Manipulate[ListPlot[QNN[[;;,l]],Joined->True,Frame->True,PlotRange->{{1,m},All},GridLines->Automatic],{l,1,mFinal,1}]

Does the DNN linear transformation contain the dominant modes, but with subdominant modes not yet eliminated?

Qred=Q[[;;,1;;mFinal]]; QredT=Transpose[Qred]; Pdom =  Qred.QredT;
Pred = IdentityMatrix[m] - Pdom;
Dimensions[Pred]

res= Simplify[Transpose[Pred.QNN]];
Norm[#]& /@ res

dom= Simplify[Transpose[Pdom.QNN]];
Norm[#]& /@ dom

DOF reduction at each stage, after an initial coordinate transformation

Define the number of degree of freedom reduction per layer δDOF,  final number of DOFs, thus defining the number of layers. Add an initial coordinate transformation layer.

lmt = NetChain[Flatten[{LinearLayer[m ,"Input"-> m ,"Biases"->None],
Table[LinearLayer[m - k δDOF,"Input"-> m - (k-1) δDOF,"Biases"->None],{k,1,nLayers}] },1]
]

Carry out SGD to find DNN weights

LinearModelt = NetTrain[lmt,data]

Define a function to extract the layer weights

Wt[l_]:=Normal[NetExtract[NetExtract[LinearModelt,l],"Weights"]];
{Dimensions[W[1]],Dimensions[W[nLayers]]}

NNtT=Wt[nLayers+1];
For[k=nLayers,k≥1,k--,NNtT=NNtT . Wt[k]];
NNt=Transpose[NNtT];
Dimensions[NNt]

Does the trained DNN carry out projection onto the SVD dominant modes? First, visualize the DNN modes

Manipulate[ListPlot[NNt[[;;,l]],Joined->True,Frame->True,PlotRange->{{1,m},All},GridLines->Automatic],{l,1,mFinal,1}]

No, from the above visualization the eigenmodes are not obtained. Orthonormalize the DNN linear transformation to check if dominant modes are obtained

{QNNtT,RNNt}=QRDecomposition[NNt]; QNNt=Transpose[QNNtT]; Dimensions[QNNt]

Manipulate[ListPlot[QNNt[[;;,l]],Joined->True,Frame->True,PlotRange->{{1,m},All},GridLines->Automatic],{l,1,mFinal,1}]

GraphicsGrid[ ArrayReshape[
Table[
ListPlot[QNNt[[;;,l]],Joined->True,Frame->True,PlotRange->{{1,m},All},GridLines->Automatic,ImageSize->Large,PlotLabel->"Mode " <> ToString[l]],{l,6}],{2,3}]]

GraphicsGrid[ ArrayReshape[
Table[
ListPlot[Q[[;;,l]],Joined->True,Frame->True,PlotRange->{{1,m},All},GridLines->Automatic,ImageSize->Large,PlotLabel->"Mode " <> ToString[l]],{l,6}],{2,3}]]

Does the DNN linear transformation contain the dominant modes, but with subdominant modes not yet eliminated?

Qred=Q[[;;,1;;mFinal]]; QredT=Transpose[Qred]; Pdom =  Qred.QredT;
Pred = IdentityMatrix[m] - Pdom;
Dimensions[Pred]

res= Simplify[Transpose[Pred.QNNt]];
Norm[#]& /@ res

dom= Simplify[Transpose[Pdom.QNNt]];
Norm[#]& /@ dom

DOF reduction at each stage, after an initial coordinate transformation, with initialization of transformation to exact solution

Define the number of degree of freedom reduction per layer δDOF,  final number of DOFs, thus defining the number of layers. Add an initial coordinate transformation layer.

QT=Transpose[Q];

lmti = NetChain[Flatten[{LinearLayer[m ,"Input"-> m ,"Biases"->None, "Weights"->QT ],
Table[LinearLayer[m - k δDOF,"Input"-> m - (k-1) δDOF,"Biases"->None],{k,1,nLayers}] },1]
]

Carry out SGD to find DNN weights

LinearModelti = NetTrain[lmti,data]

Define a function to extract the layer weights

Wti[l_]:=Normal[NetExtract[NetExtract[LinearModelti,l],"Weights"]];
{Dimensions[Wti[1]],Dimensions[Wti[nLayers]]}

NNtiT=Wti[nLayers+1];
For[k=nLayers,k≥1,k--,NNtiT=NNtiT . Wti[k]];
NNti=Transpose[NNtiT];
Dimensions[NNti]

Does the trained DNN carry out projection onto the SVD dominant modes? First, visualize the DNN modes

Manipulate[ListPlot[NNti[[;;,l]],Joined->True,Frame->True,PlotRange->{{1,m},All},GridLines->Automatic],{l,1,mFinal,1}]

No, from the above visualization the eigenmodes are not obtained. Orthonormalize the DNN linear transformation to check if dominant modes are obtained

{QNNtiT,RNNti}=QRDecomposition[NNti]; QNNti=Transpose[QNNtiT]; Dimensions[QNNti]

Manipulate[ListPlot[QNNti[[;;,l]],Joined->True,Frame->True,PlotRange->{{1,m},All},GridLines->Automatic],{l,1,mFinal,1}]

Does the DNN linear transformation contain the dominant modes, but with subdominant modes not yet eliminated?

res= Simplify[Transpose[Pred.QNNti]];
Norm[#]& /@ res

dom= Simplify[Transpose[Pdom.QNNti]];
Norm[#]& /@ dom

DOF reduction at each stage, after an initial coordinate transformation, with initialization of all layers to exact solution

Define the number of degree of freedom reduction per layer δDOF,  final number of DOFs, thus defining the number of layers. Initialize all layers to exact, SVD solution

Z[k_]:=MapThread[Join,{IdentityMatrix[k],Table[0,{i,k},{j,δDOF}]}]

Z[4] // MatrixForm

lmtia = NetChain[Flatten[{LinearLayer[m ,"Input"-> m ,"Biases"->None, "Weights"->QT ],
Table[LinearLayer[m - k δDOF,"Input"-> m - (k-1) δDOF,"Biases"->None, "Weights"->Z[m-k δDOF]],{k,1,nLayers}] },1]
]

Carry out SGD to find DNN weights

LinearModeltia = NetTrain[lmtia,data]

Define a function to extract the layer weights

Wtia[l_]:=Normal[NetExtract[NetExtract[LinearModelti,l],"Weights"]];
{Dimensions[Wtia[1]],Dimensions[Wtia[nLayers]]}

NNtiaT=Wtia[nLayers+1];
For[k=nLayers,k≥1,k--,NNtiaT=NNtiaT . Wtia[k]];
NNtia=Transpose[NNtiaT];
Dimensions[NNtia]

{QNNtiaT,RNNtia}=QRDecomposition[NNtia]; QNNtia=Transpose[QNNtiaT]; Dimensions[QNNtia]

The SGD only mildly perturbs the exact solution

Manipulate[ListPlot[{Q[[;;,l]],QNNt[[;;,l]],QNNti[[;;,l]],QNNtia[[;;,l]]},Joined->True,Frame->True,PlotRange->{{1,m},All},GridLines->Automatic,PlotLegends->{"Exact","Random","First layer","All layers"},PlotLabel->"DNN DOF reduction: initialization comparison"],{l,1,mFinal,1}]

Does the DNN linear transformation contain the dominant modes, but with subdominant modes not yet eliminated?

res= Simplify[Transpose[Pred.QNNtia]];
Norm[#]& /@ res

dom= Simplify[Transpose[Pdom.QNNtia]];
Norm[#]& /@ dom

Compare to exact

res= Simplify[Transpose[Pred.Q[[;;,1;;mFinal]]]];
Norm[#]& /@ res

dom= Simplify[Transpose[Pdom.Q[[;;,1;;mFinal]]]];
Norm[#]& /@ dom

Single linear layer

Define the number of degree of freedom reduction per layer δDOF,  final number of DOFs, thus defining the number of layers

δDOF=2; mFinal=6; nLayers=(m-mFinal)/δDOF;
l1 = NetChain[{LinearLayer[mFinal,"Input"-> m ,"Biases"->None]}]

Carry out SGD to find DNN weights

LinearModel1 = NetTrain[l1,data]

W1=Transpose[Normal[NetExtract[NetExtract[LinearModel1,1],"Weights"]]]; Dimensions[W1]

Does the trained DNN carry out projection onto the SVD dominant modes? First, visualize the DNN modes

Manipulate[ListPlot[W1[[;;,l]],Joined->True,Frame->True,PlotRange->{{1,m},All},GridLines->Automatic],{l,1,mFinal,1}]

Again, poor separation of dominant modes.

Lanczos algorithm approach

Form correlation matrix

uCT = Transpose[uC]; Corr = uC.uCT; Dimensions[Corr]

Λ=Eigenvalues[Corr]

Lanczos[n_,u_,A_]:=Module[{m,k,β0,β1,p0,p1,p2,b,v,α,T},
m=Length[A]; T = Table[0,n,n]; β0=0; 
p0=Table[0,m]; b=u; p1=b/Norm[b];
For[k=1,k≤n,k++,
v = A.p1; α = p1.v; v  = v-β0 p0-α p1;
β1 = Norm[v]; p2 = v/β1;
If[k>1,T[[k,k-1]]= β0];
T[[k,k]]=α;
If[k<n,T[[k,k+1]]=β1];
β0=β1; p0=p1; p1=p2;
];
Return[T];
]

Lanczos[5,u[[;;,1]],Corr] // MatrixForm

CompareEigvals[n_,j_]:= {Eigenvalues[Lanczos[n,u[[;;,j]],Corr]],Λ[[1;;n]]}

CompareEigvals[3,1]

CompareEigvals[5,1]

Interspersed non-linear layers