The idea here is to code the Steinhaug algorithm Alg7.2 p171 Nocedal and investigate if a fixed resolution wave form is competitive.

Steinhaug

Wave Form

I have a model function m[p]  and I want to solve the ODE
	p'[τ]=-α[τ] Dm[p[τ]]
with initial condition
	p[0]=(OverVector[0])
until I pass through the trust region constraint boundary
	||p[(SuperStar[τ])]||=Δ  or better yet p[(SuperStar[τ])].S.p[(SuperStar[τ])]==1
I want to make the progress uniform so I want to make (SuperStar[τ])=1 and
     	p'[τ].S.p[τ]==2 τ
So I want to make
	p'[τ].S.p[τ]=-α[τ] Dm[p[τ]].S.p[τ] =2τ
which makes α[τ]=-((2τ )/(Dm[p[τ]].S.p[τ]))
and I want to solve the damm thing on the interval [0,1].
to make 
	(d/(dτ))( p[τ].S.p[τ] )=(d/(dτ))[(τ^2)]
in which case p[τ].S.p[τ]=0+(τ^2) and we know when we cross the constraint boundary.  The 
auxiliary ODE is 
	p'[τ].S[τ].p[τ]=2τ
or 
	-α[τ] Dm[p[τ]].S.p[τ]=2 τ
If possible I would like to do non-spherical constraints. wave form relaxation says you start with the straight line and iterate according to
	(Subscript[p, n+1])[t]=(OverVector[0])+(∫_0^t)Dm[(Subscript[p, n])[τ]dτ
In general it needs to be implemented in quadrature with a specific fixed number of points.  For clarity and simplicity I am going to implement it with a simple quadrature hand rule and an initially uniform point distribution.

Optional Extra:
With the constraint
	1/2p[τ].S.p[τ]=(τ^2)
which differentiates to 
	p'[τ].S.p[τ] = 2 τ
I need to map stuff from [0,Δ]  to like a fixed [0,1].

Break

Clear[p,n]
F[p_]:=Module[
{dm,pp=0},
dm=Map[Dm, p];
Join[{{0,0}},
Table[
pp-= ((dm[[i]])/\[Placeholder]),(* (1/2)(dm[[i]]+dm[[i+1]]), *)
{i,1,Length[p]-1}]
]
];
p[n_]:=p[n]=F[p[n-1]]

B=RandomReal[{0,1},{2,2}];
B = B + Transpose[B];
g = RandomReal[{-1,1},2];
m[x_]:=(1/2)x.B.x - g.x
Dm[x_]:=B.x-g;
M = 12;
Δ=0.8M;
sd=-Δ((Dm[{0,0}])/(Norm[Dm[{0,0}]]))

p[0]=Table[i sd, {i, 0, 1, 1/15}];
ContourPlot[ m[{x1,x2}],{x1,-M, M},{x2, -M, M},
Epilog->{
{Red, Circle[{0,0},Δ], Point[{0,0}]},
Table[
{Hue[(i/24)],Line[p[i]], Map[Point,p[i]]},
{i,0,24}]
}
]

p[12]

p[0]
p[1]

sd=-Δ((Dm[{0,0}])/(Norm[Dm[{0,0}]]))
p[0]=Table[i sd, {i, 0, 1, 1/15}]
Dmp[0]=Map[Dm, p[0]]
p[1]=Table[
pp = p[0]

Toda Flow

To implement the Toda flow I need to split a possibly complex matrix.

SplitMatrix[B_]:=Module[
{BL = LowerTriangularize[B,-1]},
B-BL
]
SplitComplexMatrix[B_]:= Module[
{BL=LowerTriangularize[B,-1],π1,π2},
π1=(1/2)(BL+ Conjugate[BL])+ I Im[ DiagonalMatrix[Diagonal[B]]];
π2=B-π1;
{π1,π2}
]

The evolution 
	B'[t]==π1[B].B-B.π1[[B]
should be iso-spectral and should force force the lower triangular portion of
B to zero. It does seem to work.  For a real matrix the diagonal should become ordered.

Clear[F,B,t]
n=12;
TMax=300;
f[B_]:=B;
B0=RandomReal[{-1,1},{n,n}];
B0=LowerTriangularize[UpperTriangularize[B0,-1],1];
B0=B0 + Transpose[B0];
F[B_List]:= Module[
{π1},
π1=SplitComplexMatrix[B][[2]];
B.π1 - π1.B
]
BSol=B/.NDSolve[{
B'[t]==(10^2) F[B[t]],
B[0]==B0
},
B,
{t, 0, TMax},
AccuracyGoal->2,
PrecisionGoal->2
] [[1]]

The following shows that the evolution is isospectral and that the diagonal becomes ordered.

Needs["DifferentialEquations`InterpolatingFunctionAnatomy`"];
λ0=Eigenvalues[B0];
ListLogPlot[
Transpose[{
Flatten[InterpolatingFunctionGrid[BSol]],
Map[Norm[1-((Eigenvalues[#])/(λ0))]&,InterpolatingFunctionValuesOnGrid[BSol]]
}]
]

ts=Flatten[InterpolatingFunctionGrid[BSol]];
BDias=Map[Diagonal,InterpolatingFunctionValuesOnGrid[BSol]];
Fred=Table[
tt=ts[[i]];
ds=BDias[[i]];
Table[ {tt,ds[[j]]},{j, n}],
{i, Length[ts]}
];
ListLogLinearPlot[Transpose[Fred],
PlotRange->All]

OK so what we have learned.
The Toda flow really does try to set the order. In practice it does wuite well. The last pieces to resolve typically involve switching two entries which is kinda interesting.  I would guess that the knowledge of the final version could be used to advantage in a wave form relaxation algorithm.