function [result]=smoothing_spline(X_l, Y_l, X_t, Y_t, X_eval, Y_eval, lambda_list)
    X_l=X_l+1;
    X_t=X_t+1;
    X_eval=X_eval+1;
   % Dimension
   d = size(X_t,2);
   %lambda_list=2.^(-20:1:20);              % Lambda
   

    k=ceil(d/2)+1; %(ceil(d/2)+1:1:7);
    if d==1 || d==2
       k=2;
    end
      
   % noch ein paar Vorbereitungen:
   d_odd=false;
   if mod(d,2)==1
       d_odd=true;
   end
   
   % für Monomiale die Aufteilung vorbereiten und als Matrix abspeichern:
   %ausgabe(sprintf('Erstelle Monomial-Tabelle... %u Eintraege.',  nchoosek(d+max(param_k_list)-1,d)));
   
   if d>1
   
       matrix_monomials=allVL1(d,k-1,'<=');
   
   else
       
       matrix_monomials=(0:(k-1))';
       
   end
   
   % Testdaten in Cell-Array (notwendig für Parallelisierung)
   cell_test=mat2cell(X_t, ones(1,size(Y_t,1)));
   
   clear n;
   clear X_t M;
   
   %     Vorbereitung eindeutige X-Werte
   %
   % Doppelte Einträge finden, in den entsprechenden Y die Summe der Y
   % eintragen. Die Variable mit der Anzahl der Einträge zu einem Vektor i
   % heißt Y_count und ist Spaltenvektor.
   % Achtung: Der Befehl unique sortiert Reihenfolge um (aufsteigend)!

   [X_l,~,Y_l(:,2)] = unique(X_l,'rows');

   % Y-Werte zusammenrechnen:
   [N,~]=size(X_l);   %muss neu errechnet werden!
   Y_count=zeros(N,1);
   Ytmp=zeros(N,1);
   for i=1:size(Y_l)
       Ytmp(Y_l(i,2),1)=Ytmp(Y_l(i,2),1)+Y_l(i,1);
       Y_count(Y_l(i,2))=Y_count(Y_l(i,2))+1;
   end
   Y_l=Ytmp(:,1);
   clear Ytmp;
   clear i;
   
       
       % definiere l (so viele Monomiale gibt es)
        l=nchoosek(d+k-1 , d);
       
       % Theta bestimmen:
       if d_odd
           theta=gamma(d/2-k)/(2^(2*k)*pi^(d/2)*factorial(k-1));
       else
           theta=(-1)^(k+d/2+1)/(2^(2*k-1)*pi^(d/2)*factorial(k-1)*factorial(k-d/2));
       end
       
       
       
       % LGS erstellen parallel:
       % hier die Matrix-Zeilen erst als Zellen in einem Array definieren
       % ... Danach kann man sie sich wieder als Matrix konvertieren.
       % Zeiteinsparung bei einem 8-Kern-Rechner: 2/3
       
       cell_lgs=cell(N,1);
       parfor i=1:N
           
           % Eintraege der Phi-Auswertungen
           v=zeros(1,N+l);
           for j=1:l
               v(1,j+N)=prod(X_l(i,:).^matrix_monomials(j,:));
           end
       
           % Eintraege der K-Funktion:
           for j=1:N
               nz=norm(X_l(i,:)-X_l(j,:));
               if nz==0
                   v(1,j)=0;
               else
                   if d_odd
                       v(1,j)=theta*nz^(2*k-d);
                   else
                       v(1,j)=theta*nz^(2*k-d)*log(nz);
                   end
               end
           end
           
           % Erstellten Vektor als Zeile speichern
           cell_lgs{i}=v;
       end
       
       % Cell-Array umwandeln in eine Matrix
       matrix_lgs=zeros(N+l,N+l);
       matrix_lgs(1:N,1:N+l)=cell2mat(cell_lgs);
       matrix_lgs(N+1:N+l,1:N)=matrix_lgs(1:N,N+1:N+l)';
       
       % und die Anzahl der Y-Werte für einen X-Wert einbauen
       for i=1:N
           Y_count(i);
           matrix_lgs(i,1:N+l)=Y_count(i)*matrix_lgs(i,1:N+l);
       end
       
       matrix_lgs(1:N,:)=1/N*matrix_lgs(1:N,:);
       
       % Cell-Array wieder löschen
       clear cell_lgs
       %toc
         
         %ausgabe(sprintf ('Creating sparse matrix'))

         tol_factor=20*(N+l)*eps;
         spaltennorm_sq_mat=spaltennorm(matrix_lgs).^2;
         spaltennorm_sq_mat_tmp(1:N)=spaltennorm_sq_mat(1:N)+min(lambda_list)^2;
         tol=tol_factor*sqrt(max(spaltennorm_sq_mat_tmp));
         
         matrix_lgs(:,abs(spaltennorm(matrix_lgs))<tol)=0;

         % Sparse-Matrix erstellen
         [mi,mj,ms]=find(matrix_lgs);
         clear matrix_lgs;
         matrix_lgs=sparse(mi,mj,ms,N+l,N+l,size(ms,1)+N);
         clear mi mj ms;
         %toc

         %ausgabe(sprintf('Rechne Lambda durch'))
         % b-Vektor setzen: erste Einträge Y-Werte, dann Nullen
         b=1/N*Y_l;
         b(N+1:N+l,1)=0;
         
         % Lambdas durchrechnen:
         number_lambdas=length(lambda_list');

         % DMatrix f. Lamda addieren
         mat_diag=sparse(1:N,1:N,ones(1,N),N+l,N+l);
         
         %opts.solution='min2norm';
         %ausgabe (sprintf('Starte Rechnungen ...'))
         parfor i=1:number_lambdas
             %tic
             % Lambda einsetzen (geht so am schnellsten)
             matrix_lgs_tmp=matrix_lgs+lambda_list(i)*mat_diag;
   
             % das LGS loesen und als Zeile in die Matrix B schreiben:
             spaltennorm_sq_mat_tmp=spaltennorm_sq_mat(1:N)+lambda_list(i)^2;
             tol=tol_factor*sqrt(max(spaltennorm_sq_mat_tmp));
             warning('off','MATLAB:rankDeficientMatrix');
             B(i,:)=spqr_solve(matrix_lgs_tmp , b ,struct ('solution','min2norm','tol',tol));%B(i,:)=spqr_solve(matrix_lgs_tmp , b ,opts,'tol',tol);  %,'tol',10^(-5)))';
             warning('on','MATLAB:rankDeficientMatrix');
             %ausgabe (sprintf('k=%u, L=%u,  time needed: %u', param_k_list(i_loop_k), lambda_list(i), %toc ))
         end

         % das LGS brauchen wir nicht mehr ... hoffe ich ;)      
         clear matrix_lgs b;
            
          vec_error=zeros(size(B,1),1);
          parfor i=1:size(Y_t,1) 
              p=zeros(N+l,1);
              for l_count=1:l
                  p(N+l_count)=prod(cell_test{i}.^matrix_monomials(l_count,:));
              end
              for N_count=1:N
                  nz=norm(cell_test{i}-X_l(N_count,:),2);
                  if nz==0
                      p(N_count)=0;
                  else
                      if d_odd
                          p(N_count)=theta*nz^(2*k-d);
                      else
                          p(N_count)=theta*nz^(2*k-d)*log(nz);
                      end
                  end
              end
    
              vec_error=vec_error+(max(B*p,0)-Y_t(i,1)).^2;
          end
          %ausgabe('Werte berechnet')
         
         
         % ein paar Daten
         [~,best_lambda]=min(vec_error);
         bestB=B(best_lambda,:);
         
         clear B;
         %ausgabe(sprintf ('Best error for k=%u: %u with Lambda=%u',  param_k_list(i_loop_k), best_error(i_loop_k,1), lambda_list(best_error(i_loop_k,2))))
   
   
%    clear l theta;
%    [~,best_err_line] = min(best_error(:,1));
%    %save(filename, 'best_err_line', '-append');
%    
%    
%    %ausgabe(sprintf('Bester Fehler gefunden: %u bei k=%u und Lambda=%u', best_error(best_err_line,1), param_k_list(best_err_line), lambda_list(best_err_line)))
% %   sqrt(best_error(best_err_line,1))/N/param_samplerate
%    
%    bestB=bestB{best_err_line};
   
   
   
   
   % Daten durchrechnen:
   %ausgabe('Evaluationsdaten durchrechnen');
   %load(filename, 'X_eval', 'Y_eval')
%    k=param_k_list(best_err_line);
%    l=nchoosek(d+k-1 , d);
   
%    % Theta
%    if d_odd
%        theta=gamma(d/2-k)/(2^(2*k)*pi^(d/2)*factorial(k-1));
%    else
%        theta=(-1)^(k+d/2+1)/(2^(2*k-1)*pi^(d/2)*factorial(k-1)*factorial(k-d/2));
%    end
   
   cell_eval=mat2cell(X_eval, ones(1,size(X_eval,1)));
   
   
   parfor i=1:size(X_eval,1)
       p=zeros(N+l,1);
       for l_count=1:l
           p(N+l_count)=prod(cell_eval{i}.^matrix_monomials(l_count,:));
       end
       for N_count=1:N
           nz=norm(cell_eval{i}-X_l(N_count,:),2);
           if nz==0
               p(N_count)=0;
           else
               if d_odd
                   p(N_count)=theta*nz^(2*k-d);
               else
                   p(N_count)=theta*nz^(2*k-d)*log(nz);
               end
           end
       end
       Y_calc(i,1)=max(bestB*p,0);
   end
   
    result.mse=mean((Y_calc-Y_eval).^2);
    result.res=Y_calc-Y_eval;
    result.param=lambda_list(best_lambda);

end




function v=spaltennorm(X)
   v=zeros(1,size(X,2));
   for i=1:size(X,2)
       v(i)=norm(X(:,i),2);
   end
end