n-Array Inverses
Full code
n-Array Inverses
Full code
Full code
I. Basic Symbols
I.0. Memory Cleanup
I.0. Memory Cleanup
In[]:=
"Evaluation Date:"DateObject[]Remove["Global`*"]
Out[]=
Evaluation Date:
Out[]=
I.1. Symbols & Patterns
I.1. Symbols & Patterns
I.1 Symbol definitions
I.1 Symbol definitions
In[]:=
"1) Assumptions and definitions""Evaluation Date:"DateObject[]"1.0) Symbol patterns:";head[h_,c_]:=component[H_,di_]:=tensor[H_,dimsA_]:="1.1) Assumptions patterns:";assDim[dimsA_,di_]:=FullSimplify[Element[dimsA|di,PositiveIntegers]&&(Thread[di<dimsA]/.List->And)]assT[H_,h_,dimsA_]:=Element[tensor[head[H,h],dimsA],Arrays[dimsA,Reals]]assC[H_,h_,di_]:=Element[component[head[H,h],di],NonNegativeReals]"2.0) Symbol examples:"Dk={D1,D2};ik={i1,i2};"head"->head[H,h]"dimensions"->Dk&&ik"tensor component"->component[head[H,h],ik]"abstract tensor"->tensor[head[H,h],Dk]"2.1) Assumptions examples:""dimension assumptions"->assDim[Dk,ik]"tensor assumptions"->assT[H,h,Dk]"component assumptions"->assC[H,h,ik]
h
c
H
di
H
dimsA
Out[]=
1) Assumptions and definitions
Out[]=
Evaluation Date:
Out[]=
Out[]=
2.0) Symbol examples:
Out[]=
head
H
h
Out[]=
dimensions{D1,D2}&&{i1,i2}
Out[]=
tensor component
H
h
{i1,i2}
Out[]=
abstract tensor
H
h
{D1,D2}
Out[]=
2.1) Assumptions examples:
Out[]=
dimension assumptions(D1|D2|i1|i2)∈&&i1>0&&i2>0&&i1<D1&&i2<D2
Out[]=
tensor assumptions∈Arrays[{D1,D2},,{}]
H
h
{D1,D2}
Out[]=
component assumptions∈&&≥0
H
h
{i1,i2}
H
h
{i1,i2}
II.1 Calculation Patterns: array inverses
II.1 Calculation Patterns: array inverses
REFERENCES TO ELABORATE THE ALGORITHM:
1) M. Liang and B. Zheng ; “Further results on Moore–Penrose inverses of tensors”; Computers and Mathematics with Applications 77 (2019) 1282–1293
2) Brett W. Bader and Tamara G. Kolda ; “MATLAB Tensor Classes for Fast Algorithm Prototyping”; SANDIA REPORT SAND2004-5187 (October 2004)
II.1.1 Indexing arrays
II.1.1 Indexing arrays
In[]:=
"Evaluation Date:"DateObject[]"A) Lists of indexes";SubSetsLength[ss1_,ss2_]:=Map[Length,{ss1,ss2}]SubSetsDimensions[dimA_,ss1_,ss2_]:={dimA[[ss1]],dimA[[ss2]]}ArrayPositions[dimA_]:=Module[{dA,dAi,range},dA=Length[dimA];dAi=Table[d[i],{i,dA}];range=Table[{d[i],dimA[[i]]},{i,dA}];Flatten[Apply[Table,Flatten[{{{dAi}},range},1]],dA]]PositionSubSets[dimA_,ss1_,ss2_]:=Table[{ArrayPositions[dimA][[iN,ss1]],ArrayPositions[dimA][[iN,ss2]]},{iN,Length[ArrayPositions[dimA]]}]MatrixIndex[dimA_,ss1_,ss2_,iN_]:=Module[{dK,dL,dimK,dimL,iK,iL,j1,j2},dK=Length[ss1];dL=Length[ss2];dimK=SubSetsDimensions[dimA,ss1,ss2][[1]];dimL=SubSetsDimensions[dimA,ss1,ss2][[2]];iK=PositionSubSets[dimA,ss1,ss2][[iN,1]];iL=PositionSubSets[dimA,ss1,ss2][[iN,2]];j1=iK[[1]]+Sum[(iK[[dk]]-1)*Product[dimK[[jK]],{jK,dk-1}],{dk,2,dK}];j2=iL[[1]]+Sum[(iL[[dl]]-1)*Product[dimL[[jL]],{jL,dl-1}],{dl,2,dL}];{j1,j2}]
Out[]=
Evaluation Date:
Out[]=
II.1.2 Folding and Unfolding
II.1.2 Folding and Unfolding
In[]:=
"Evaluation Date:"DateObject[]"A) UnFolding maps";IndexMap[dimA_,ss1_,ss2_,iN_]:={ArrayPositions[dimA][[iN]],MatrixIndex[dimA,ss1,ss2,iN]}UnfoldArray[A_,dimA_,ss1_,ss2_]:=SparseArray[Table[IndexMap[dimA,ss1,ss2,iN][[2]]->Extract[A,IndexMap[dimA,ss1,ss2,iN][[1]]],{iN,Apply[Times,dimA]}]]//Normal"B) Reverse UnFolding maps";ReverseIndexMap[dimA_,ss1_,ss2_,iN_]:={MatrixIndex[dimA,ss1,ss2,iN],ArrayPositions[dimA][[iN]]}FoldArray[M_,dimA_,dimM_,ss1_,ss2_]:=SparseArray[Table[ReverseIndexMap[dimA,ss1,ss2,iN][[2]]->Extract[M,ReverseIndexMap[dimA,ss1,ss2,iN][[1]]],{iN,Apply[Times,dimM]}]]//Normal
Out[]=
Evaluation Date:
Out[]=
II.1.3 Array inverses
II.1.3 Array inverses
In[]:=
"Evaluation Date:"DateObject[]InvertedSubSets[ss1_,ss2_]:=TakeList[Table[i,{i,Length[ss1]+Length[ss2]}],{Length[ss2],All}];unfoldedPseudoInverse[A_,ss1_,ss2_]:=Module[{ass,dimA,m},ass=Element[DeleteCases[Flatten[A],0],PositiveReals];dimA=Dimensions[A];m=UnfoldArray[A,dimA,ss1,ss2];Assuming[ass,Refine[PseudoInverse[m]]]]foldedPseudoInverse[A_,ss1_,ss2_]:=Module[{uPINV,dimA,dimuPINV,perm,dimPINV,iss},uPINV=unfoldedPseudoInverse[A,ss1,ss2];dimA=Dimensions[A];dimuPINV=Dimensions[uPINV];dimPINV=Flatten[{dimA[[ss2]],dimA[[ss1]]}];iss=InvertedSubSets[ss1,ss2];FoldArray[uPINV,dimPINV,dimuPINV,iss[[1]],iss[[2]]]]
Out[]=
Evaluation Date:
Out[]=
II.1.4 Array inner product
II.1.4 Array inner product
"Evaluation Date:"DateObject[]ArrayInner[A1_,A2_,slots_]:=Inactive[FullSimplify][Inactive[TensorContract][Inactive[TensorProduct][A1,A2],slots]]slotsAXA1[ss2_]:=Table[{ss2[[i]],Length[ss2]+ss2[[i]]},{i,Length[ss2]}]slotsAXA2[ss1_]:=Table[{InvertedSubSets[ss1,ss1][[2,i]],Length[ss1]+InvertedSubSets[ss1,ss1][[2,i]]},{i,Length[InvertedSubSets[ss1,ss1][[2]]]}]slotsXAX1[ss1_,ss2_]:=Table[{Length[ss2]+ss1[[i]],Length[ss1]+Length[ss2]+ss1[[i]]},{i,Length[ss1]}]slotsXAX2[ss1_,ss2_]:=Table[{InvertedSubSets[ss2,ss2][[2,i]],Length[ss2]+InvertedSubSets[ss2,ss2][[2,i]]},{i,Length[InvertedSubSets[ss2,ss2][[2]]]}]
II. Testing the algorithm
II.1) Array Definitions:
II.1) Array Definitions:
II.2) Tests:
II.2) Tests:
II.2.1) Folding indexes:
II.2.1) Folding indexes:
II.2.2) Fold/Unfold:
II.2.2) Fold/Unfold:
II.2.3) Inverse arrays:
II.2.3) Inverse arrays:
II.2.3) Testing if those are really generalized inverses:
II.2.3) Testing if those are really generalized inverses: