kfli

30 Reputation

4 Badges

3 years, 233 days

MaplePrimes Activity


These are replies submitted by kfli

@Carl Love 

Yes, that works. I'm using Maple 2018. Any reason for the compound line other than aesthetic? 

I must have made a mistake before because I ran the different alternatives again and you are correct, Alt1 and Alt3 are the same. To be honest I have no idea why it was showing different result before. It was a simple copy paste and print result.

@kfli 
This seems to be correct also

R := norm(df,2);
Q := df/~R;
gamma := Transpose(Q).fval/~R;
 

@Joe Riel

That's really great! I will be using that in all my codes from now on :D I was really missing the Matlab profiling function in Maple.  

@Carl Love 

Hah that's why I need people like you. I was ripping my hair out over that and it wasn't even the issue. 

Well the code does 'work' in the sense that it gives the right solution to the set of equations. There was only that weird issue where the produre just stopped for no apparant reason. No errors or anything, and I thought it had to do with the procedure issue but I was mistaken. I changed a few things as you see. I changed QR deomposition to 'NAG' output and I check the condition number of Q instead of R (because R is now a vector). I am not sure if that is the same thing mathematically so it could be working slightly different than original, but it atleast finishes the procedure. 

I might have been confused on the line of the code you comment. I was trying to implement MTM[ldivide]. I tried all three alternatives.

Alt.1 
R := norm(df,2);
Q := R/~df;
gamma := R/~Transpose(Q).fval;


Alt. 2
R := norm(df,2);
Q := MTM[ldivide](R,df);
gamma := MTM[ldivide](R,Transpose(Q).fval);


Alt. 3
gamma := Transpose(df).fval;

All give different results (ldivide is corrent). I thought that Alt.1 and Alt.2 were equivalent but I was wrong. Alt.2 is the from the original Matlab code. I thought lvdivide was the same as '/~'. The idea is to solve df.gamma=fval with QR decomp. 

df.gamma=fval
QR.gamma=fval
gamma=R^-1.Q' . fval

Side note: If we can make a procedure that alters prexisting Q R instead of making QR from scratch everytime we change df than this code would be working smooth as eggs.

@Carl Love 
That's crazy. Thanks for the info. Where did you read this? As of right now I have changed some things so that the code works for all types of equations. It looks like this if you are interested.

 

restart:
Digits:=12:
Runtime:=time():
with(LinearAlgebra):
with(plots):
with(orthopoly):
with(ArrayTools):
    


N := 2:
phix := Vector(N):
X    := Vector(N):
phix[1] := cos(x[2]):
phix[2] := 3*cos(x[1]):
X[1] := 1.0:
X[2] := 2.0:

#Code AndersonAcceleration.
AndersonAcceleration:=proc(N,phi,X)
global x, xS, here1, here2;
local mMax, itmax, atol, rtol, droptol, beta, AAstart, res_hist, df, DGg, gamma, X0;
local DG, mAA, iter, gval,fval,res_norm, condDF, tol, f_old, g_old, y, i, k, Q, R, QRg, dfT, DGT;

(*
%------------------------------------------------------------------------
% Function xS = AndersonAcceleration(N,phi,x0).
%
% Fixed-point iteration with or without Anderson acceleration.
% 'phi' is the fixed point iteration map and 'xS' is the 
% solution, so that xS = phi(xS).
%
% Input arguments:
%    X0 = Initial value solution. Note: In this function this variable 
%        must be a column vector.
%       1. 'mMax' = maximum number of stored residuals (non-negative integer).
%       NOTE: 'mMax' = 0 => no acceleration. default=1000
%       2. 'itmax' = maximum allowable number of iterations. default=1000
%       3. 'atol' = absolute error tolerance. default=1.0e-6
%       4. 'rtol' = relative error tolerance. default=1.0e-3
%       5. 'droptol' = tolerance for dropping stored residual vectors to 
%       improve conditioning: If 'droptol' > 0, drop residuals if the
%       condition number exceeds droptol; if droptol <= 0,
%       do not drop residuals.
%       6. 'beta' = damping factor: If 'beta' > 0 (and beta ~= 1), then 
%       the step is damped by beta; otherwise, the step is not damped.
%       NOTE: 'beta' can be a function handle; form beta(iter), where iter 
%       is the iteration number and 0 < beta(iter) <= 1.
%       7. 'AAstart' = acceleration delay factor: If 'AAstart' > 0, start 
%       acceleration when iter = AAstart.
%
% Output:
% xS = Solution vector.
%
% The notation used is that of H.F. Walker: Anderson Acceleration:
% Algorithms and implementation
%------------------------------------------------------------------------
*)

mMax    := 100:
itmax   := 100:
atol    := 1.0e-8:
rtol    := 1.0e-6:
droptol := 1.0e12:
beta    := 1.0:
AAstart := 0:

# Initialize storage arrays and number of stored residuals.
DG := Matrix():
df := Matrix():
DGg := Vector(N);
QRg := Vector(N);
mAA := 0:

X0 := X;
for iter from 0 to itmax do

   x:=X0:
   gval := Vector(phi):
   fval := gval - X0:
   res_norm := norm(fval,2):

   # Set the residual tolerance on the initial iteration.
   if iter = 0 then
      tol := max(atol,rtol*res_norm):
   fi:
    
   # Convergence test, if converged the loop stops.
   if res_norm <= tol then
      print('iter',iter,res_norm);
      break;   # Breaks for-loop
   fi:
    
   # If resnorm is larger than 1e8 at iter > 5, problem stops
   if res_norm >1e8 and iter > 5 then
      print('no_convergence',res_norm);
      break; # Breaks for-loop, diverged
   fi:

   # Fixed point iteration without acceleration, if mMax == 0.
   if mMax = 0 or iter < AAstart then
      # We update E <- g(E) to obtain the next approximate solution.
      X0 := gval:     
   else
      # With Anderson acceleration.
      # Update the df vector and the DG array.
      if iter > AAstart then
         if mAA < mMax or Size(df,2) = 1 then
            df := Concatenate(2,df,fval-f_old):
            DG := Concatenate(2,DG,gval-g_old):
         else 
            df := Concatenate(2,df[..,-1],fval-f_old):
            DG := Concatenate(2,DG[..,-1],gval-g_old):   
         fi:
         mAA := mAA + 1:
      fi:   # iter
      # We define the old g and f values for the next iteration
      f_old := fval;
      g_old := gval;

      if mAA = 0 then
         # Initialization
         # If mAA == 0, update X <- g(X) to obtain themah next approximate
         # solution. No least-squares problem is solved for iter = 0
         X0 := gval;
      else
         if mAA > 1 then
            Q,R := QRDecomposition(df,datatype=float,output='NAG');
            if type(Q, 'Matrix'(square)) then
              condDF := ConditionNumber(Q);
            else
              condDF := 1;
            fi:  
            while condDF > droptol and mAA > 1 do
                if mAA = 2 then
                   df := convert(df[..,2..-1],Vector);
                   DG := convert(DG[..,2..-1],Vector);
                else
                   df := df[..,2..-1];
                   DG := DG[..,2..-1];
                fi:
                Q,R := QRDecomposition(df,datatype=float,output='NAG');
                mAA := mAA - 1;
                if type(Q, 'Matrix'(square)) then
                  condDF := ConditionNumber(Q);
                else
                  condDF := 1;
                fi:
            od:
            if Size(df,2) > 1 then
               gamma := LeastSquares([Q,R],fval);
            else
               R := norm(df,2);
               Q := R/~df;
               gamma := R/~Transpose(Q).fval;
            fi:
            
         else
            R := norm(df,2);
            Q := R/~df;
            gamma := R/~Transpose(Q).fval;
         fi:
 
         if Size(gamma,1) > 1 then
            DGg:=DG.gamma:
         else
            DGg:=DG*gamma;
         fi:
         
         # Update the approximate solution.
         X0 := gval - DGg;

         # Damping for non-zero beta
         if beta > 0 and beta <> 1 then
            if mAA = 1 then
               QRg := Q*R*gamma;
            else
               QRg := df.gamma;
            fi:
            X0 := X0 - (1-beta)*(fval - QRg);
         fi:# isa(beta ...
      fi: # mAA = 0
   fi:# mMax == 0

od:
X[1..N] := X0[1..N];
end:

AndersonAcceleration(N,phix,X):
 

@tomleslie 

Hey thank you for the tips! The specific example above did not work on my computer. For some reason it just stopped. No error, no if statement break or nothing. So I was hoping for a solution to that but if it worked for you then I am not sure what is wrong. 

#1 Yes you are right, that was a bit lazy I suppose. I will fix that. 
I used Q := R/~df;

#2 I can't believe I didn't think of that... since I had the same notations in other parts of the code.

#3 Thanks, that's  fixed now.

@Adri van der Meer 

Right, is there a way to check in an if-statement if R1 is

R1 := [2];
 or 
R1 := 2;

for example?

And is there a way to convert R1 into a regular scalar instead of a 1x1 matrix?

@tomleslie 

 

Hey! 

Yes that was what I was looking for... but 

My specific procedure needs to be able to have differing x vector length

for example 

restart;

x = [x1 x2... xN]
myV:=proc(x)

                    N=length(x);
                    Vector(N,[x1+2*x2, x2^2+x1,..., xN+x[1]]);
          end proc;
CodeGeneration:-Matlab( myV,
                                            output="C:/Users/TomLeslie/Desktop/test.m"
                                         );

 

something like this... would that work you think?

@Carl Love @mehdi jafari 

Thanks for the reply!
I need the matrix for evaulations in another code. I build the equations and Jacobian in one file, save them, then upload them in another file wheere they are repeatedly evaluated/solved.

The jacobian looks something like 

J[1,1] = B*(x[1]+x[2])
J[1,2] = B*(3*x[2]) ....
etc. This is just made up but that is basically what the elements look like.

Then in the "solver" file I compute B and then evaluate the J with B and x values. So basically I only need it for numeric evaluations, not symbolic analysis.
 

@Carl Love Thank you for the reply! I was hoping the fix was going to be a easier. Do you know of any way I could apply the freeze and thaw to my specific problem? I couldn't seem to get it work.

Page 1 of 1