基于SMO方法的支持向量机Pascal代码实现

SMO_SVM.pas：

{
Author: Liu Liu
Website: www.aivisoft.net
E-Mail: [e-mail]
This code is translated from a p-code by J. Platt.
This code shows a SMO(Sequential Minimal Optimization) SVM.
Enjoy fast and efficient SVM!
This code now provide 4 different kernels.
}
unit SMO_SVM;

interface

uses
  Windows, SysUtils, Variants, Classes, Graphics, Math, unitypes;

type
  TVector = record
    Items: TSingleExtendedArray;
    Target: Extended;
  end;

  TVectorArray = array of TVector;

  TKernels = (klLinear, klGsRBF, klExpRBF, klPolynomial, klTanh);
  {
  klLinear: Linear Kernel
  klGsRBF: Gaussian RBF Kernel
  klExpRBF: Exponential RBF Kernel
  klPolynomial: Polynomial Kernel
  klTanh: Tanh Kernel
  }
  TSVM = class
  private
    n_Degree, N: longint;
    tolerance, eps, b: Extended;
    alph, Error_Cache, w, Precomputed_Self_Dot_Product: TSingleExtendedArray;
    function Dot_Compute(i1, i2: longint): Extended;
    function Learned_Func(k: longint): Extended;
    function Kernel_Func(i1, i2: longint): Extended;
    function PKernel_Func(i1: longint; Items2: TSingleExtendedArray): Extended;
    function ExamineExamples(i2: longint): longint;
    function TakeStep(i1, i2: longint): boolean;
  public
    end_support_i: longint; Vectors: TVectorArray;
    C, two_sigma_squared, tha, thb, pyp: Extended;
    Kernel: TKernels;
    procedure Learn;
    procedure Init(sn_Degree: longint; sC, stwo_sigma_squared, stha, sthb, spyp: Extended; sKernel: TKernels);
    procedure LearnExamples(Data: TSingleExtendedArray; Target: Extended);
    function Optimize: longint;
    function Predict(Itemk: TSingleExtendedArray): Extended;
    function SaveToFile(FileName: string): boolean;
    function LoadFromFile(FileName: string): boolean;
    destructor Destroy; override;
  end;

implementation

function TSVM.Dot_Compute(i1, i2: longint): Extended;
var
  i: longint;
  dot: Extended;
begin
  dot := 0;
  for i := 0 to n_Degree - 1 do
    dot := dot + Vectors[i1].Items[i] * Vectors[i2].Items[i];
  Result := dot;
end;

procedure TSVM.Init(sn_Degree: longint; sC, stwo_sigma_squared, stha, sthb, spyp: Extended; sKernel: TKernels);
var
  i: longint;
begin
  n_Degree := sn_Degree;
  setlength(Vectors, $100);
  C := sC; eps := 0.001; tolerance := 0.001;
  two_sigma_squared := stwo_sigma_squared;
  tha := stha; thb := sthb; pyp := spyp;
  N := 0; Kernel := sKernel;
  setlength(w, n_Degree);
end;

procedure TSVM.LearnExamples(Data: TSingleExtendedArray; Target: Extended);
begin
  if N > High(Vectors) then
    setlength(Vectors, N + $100);
  setlength(Vectors[N].Items, n_Degree);
  Vectors[N].Items := Data;
  Vectors[N].Target := Target;
  Inc(N);
end;

function TSVM.Learned_Func(k: longint): Extended;
var
  s: Extended;
  i: longint;
begin
  s := 0;
  case Kernel of
    klLinear: begin
        for i := 0 to n_Degree - 1 do
          s := s + w[i] * Vectors[k].Items[i];
      end;
  else begin
      for i := 0 to end_support_i - 1 do
          s := s + alph[i] * Vectors[i].Target * Kernel_Func(i, k);
    end;
  end;
  s := s - b;
  Result := s;
end;

function TSVM.Kernel_Func(i1, i2: longint): Extended;
var
  s: Extended;
  i: longint;
begin
  s := 0;
  case Kernel of
    klLinear: begin
        for i := 0 to n_Degree - 1 do
          s := s + Vectors[i1].Items[i] * Vectors[i2].Items[i];
        Result := s;
      end;
    klPolynomial: begin
        for i := 0 to n_Degree - 1 do
          s := s + Vectors[i1].Items[i] * Vectors[i2].Items[i];
        Result := Exp(Ln(s + 1) * pyp);
      end;
    klExpRBF: begin
        s := Precomputed_Self_Dot_Product[i1] + Precomputed_Self_Dot_Product[i2] - 2 * Dot_Compute(i1, i2);
        Result := exp(-sqrt(s) / Two_Sigma_Squared);
      end;
    klGsRBF: begin
        s := Precomputed_Self_Dot_Product[i1] + Precomputed_Self_Dot_Product[i2] - 2 * Dot_Compute(i1, i2);
        Result := exp(-s / Two_Sigma_Squared);
      end;
    klTanh: begin
        for i := 0 to n_Degree - 1 do
          s := s + Vectors[i1].Items[i] * Vectors[i2].Items[i];
        Result := Tanh(tha * s + thb);
      end;
  end;
end;

function TSVM.ExamineExamples(i2: longint): longint;
var
  y2, alph2, E2, r2: Extended;
  k, k0, i1: longint;
  tmax, E1, temp: Extended;
begin
  y2 := Vectors[i2].Target;
  alph2 := alph[i2];

  if (alph2 > 0) and (alph2 < C) then
    E2 := Error_Cache[i2]
  else
    E2 := Learned_Func(i2) - y2;

  r2 := y2 * E2;
  if ((r2 < -tolerance) and (alph2 < C)) or ((r2 > tolerance) and (alph2 > 0)) then begin
    i1 := -1; tmax := 0;
    for k := 0 to end_support_i - 1 do begin
      if (alph[k] > 0) and (alph[k] < C) then begin
        E1 := Error_Cache[k];
        temp := Abs(E1 - E2);
        if (temp > tmax) then begin
          tmax := temp;
          i1 := k;
        end;
      end;
    end;

    if i1 >= 0 then
      if takeStep(i1, i2) then begin
        Result := 1;
        exit;
      end;

    Randomize;
    k0 := Random(end_support_i);
    for k := k0 to end_support_i + k0 - 1 do begin
      i1 := k mod end_support_i;
      if (alph[i1] > 0) and (alph[i1] < C) then
        if (takeStep(i1, i2)) then begin
          Result := 1;
          exit;
        end;
    end;

    Randomize;
    k0 := Random(end_support_i);
    for k := k0 to end_support_i + k0 - 1 do begin
      i1 := k mod end_support_i;
      if takeStep(i1, i2) then begin
        Result := 1;
        exit;
      end;
    end;
  end;

  Result := 0;
end;

function TSVM.TakeStep(i1, i2: longint): boolean;
var
  i: longint;
  y1, y2, s: Extended;
  alph1, alph2, gamma, c1, c2, b1, b2, bnew, delta_b: Extended;
  a1, a2: Extended;
  E1, E2, L, H, k11, k22, k12, eta, Lobj, Hobj, t, t1, t2: Extended;
begin
  if (i1 = i2) then begin
    Result := false;
    exit;
  end;

  alph1 := alph[i1];
  y1 := Vectors[i1].Target;
  if (alph1 > 0) and (alph1 < C) then
    E1 := Error_Cache[i1]
  else
    E1 := Learned_Func(i1) - y1;

  alph2 := alph[i2];
  y2 := Vectors[i2].Target;
  if (alph2 > 0) and (alph2 < C) then
    E2 := Error_Cache[i2]
  else
    E2 := Learned_Func(i2) - y2;

  s := y1 * y2;

  if (y1 = y2) then begin
    gamma := alph1 + alph2;
    if (gamma > C) then begin
      L := gamma - C;
      H := C;
    end else begin
      L := 0;
      H := gamma;
    end;
  end else begin
    gamma := alph1 - alph2;
    if (gamma > 0) then begin
      L := 0;
      H := C - gamma;
    end else begin
      L := -gamma;
      H := C;
    end;
  end;

  if (L = H) then begin
    Result := false;
    exit;
  end;

  k11 := Kernel_Func(i1, i1);
  k12 := Kernel_Func(i1, i2);
  k22 := Kernel_Func(i2, i2);
  eta := 2 * k12 - k11 - k22;

  if (eta < 0) then begin
    a2 := alph2 + y2 * (E2 - E1) / eta;
    if (a2 < L) then a2 := L
    else if (a2 > H) then a2 := H;
  end else begin
    c1 := eta / 2;
    c2 := y2 * (E1 - E2) - eta * alph2;
    Lobj := c1 * L * L + c2 * L;
    Hobj := c1 * H * H + c2 * H;
    if (Lobj > Hobj + eps) then a2 := L
    else if (Lobj < Hobj - eps) then a2 := H
    else a2 := alph2;
  end;

  if a2 < 1E-8 then a2 := 0
  else if a2 > C - 1E-8 then a2 := C;

  if abs(a2 - alph2) < eps * (a2 + alph2 + eps) then begin
    Result := false;
    exit;
  end;

  a1 := alph1 - s * (a2 - alph2);
  if (a1 < 0) then begin
    a2 := a2 + s * a1;
    a1 := 0;
  end else if (a1 > C) then begin
    t := a1 - C;
    a2 := a2 + s * t;
    a1 := C;
  end;

  if (a1 > 0) and (a1 < C) then
    bnew := b + E1 + y1 * (a1 - alph1) * k11 + y2 * (a2 - alph2) * k12
  else begin
    if (a2 > 0) and (a2 < C) then
      bnew := b + E2 + y1 * (a1 - alph1) * k12 + y2 * (a2 - alph2) * k22
    else begin
      b1 := b + E1 + y1 * (a1 - alph1) * k11 + y2 * (a2 - alph2) * k12;
      b2 := b + E2 + y1 * (a1 - alph1) * k12 + y2 * (a2 - alph2) * k22;
      bnew := (b1 + b2) / 2;
    end;
  end;

  delta_b := bnew - b;
  b := bnew;

  t1 := y1 * (a1 - alph1);
  t2 := y2 * (a2 - alph2);

  if Kernel = klLinear then
    for i := 0 to n_Degree - 1 do
      w[i] := w[i] + t1 * Vectors[i1].Items[i] + t2 * Vectors[i2].Items[i];

  for i := 0 to end_support_i - 1 do
    if (0 < alph[i]) and (alph[i] < C) then
      Error_Cache[i] := Error_Cache[i] + t1 * Kernel_Func(i1, i) + t2 * Kernel_Func(i2, i) - delta_b;
  Error_Cache[i1] := 0;
  Error_Cache[i2] := 0;

  alph[i1] := a1;
  alph[i2] := a2;

  Result := true;
end;

procedure TSVM.Learn;
var
  numChanged, i: longint;
  examineAll: boolean;
begin
  end_support_i := N;
  setlength(alph, N); setlength(Error_Cache, N);
  setlength(Precomputed_Self_Dot_Product, N);
  for i := 0 to n_Degree - 1 do w[i] := 0;
  for i := 0 to N - 1 do begin
    alph[i] := 0;
    Precomputed_Self_Dot_Product[i] := Dot_Compute(i, i);
  end;
  b := 0;
  numChanged := 0;
  examineAll := true;
  while (numChanged > 0) or examineAll do begin
    numChanged := 0;
    if examineAll then begin
      for i := 0 to N - 1 do
        Inc(numChanged, ExamineExamples(i));
    end else begin
      for i := 0 to N - 1 do
        if ((alph[i] > 1E-8) or (alph[i] < -1E-8)) and ((alph[i] > C + 1E-8) or (alph[i] < C - 1E-8)) then
          Inc(numChanged, ExamineExamples(i));
    end;
    if examineAll then
      examineAll := false
    else if (numChanged = 0) then
      examineAll := true;
  end;
end;

function TSVM.PKernel_Func(i1: longint; Items2: TSingleExtendedArray): Extended;
var
  s: Extended;
  i: longint;
begin
  s := 0;
  case Kernel of
    klLinear: begin
        for i := 0 to n_Degree - 1 do
          s := s + Vectors[i1].Items[i] * Items2[i];
        Result := s;
      end;
    klPolynomial: begin
        for i := 0 to n_Degree - 1 do
          s := s + Vectors[i1].Items[i] * Items2[i];
        Result := Exp(Ln(s + 1) * pyp);
      end;
    klExpRBF: begin
        for i := 0 to n_Degree - 1 do
          s := s + sqr(Vectors[i1].Items[i] - Items2[i]);
        Result := exp(-sqrt(s) / Two_Sigma_Squared);
      end;
    klGsRBF: begin
        for i := 0 to n_Degree - 1 do
          s := s + sqr(Vectors[i1].Items[i] - Items2[i]);
        Result := exp(-s / Two_Sigma_Squared);
      end;
    klTanh: begin
        for i := 0 to n_Degree - 1 do
          s := s + Vectors[i1].Items[i] * Items2[i];
        Result := tanh(tha * s + thb);
      end;
  end;
end;

function TSVM.Optimize: longint;
var
  Change: boolean;
  i, j: longint;
begin
  Change := true;
  while Change do begin
    Change := false;
    for i := 0 to end_support_i - 1 do
      if (alph[i] < 1E-8) and (alph[i] > -1E-8) then begin
        for j := i + 1 to end_support_i - 1 do begin
          Vectors[j - 1] := Vectors[j];
          alph[j - 1] := alph[j];
        end;
        Change := true;
        Dec(end_support_i);
        break;
      end;
  end;
  N := end_support_i;
  setlength(Vectors, N); setlength(alph, N);
  Result := end_support_i;
end;

function TSVM.Predict(Itemk: TSingleExtendedArray): Extended;
var
  s: Extended;
  i: longint;
begin
  s := 0;
  case Kernel of
    klLinear: begin
        for i := 0 to n_Degree - 1 do
          s := s + w[i] * Itemk[i];
      end;
  else begin
      for i := 0 to end_support_i - 1 do
          s := s + alph[i] * Vectors[i].Target * PKernel_Func(i, Itemk);
    end;
  end;
  s := s - b;
  Result := s;
end;

function TSVM.SaveToFile(FileName: string): boolean;
var
  i, j: longint;
  SaveStream: TMemoryStream;
begin
  SaveStream := TMemoryStream.Create;
  SaveStream.Seek(0, soFromBeginning);
  try
    SaveStream.Write(end_support_i, sizeof(end_support_i));
    SaveStream.Write(n_Degree, sizeof(n_Degree));
    for i := 0 to end_support_i - 1 do begin
      SaveStream.Write(Vectors[i].Target, sizeof(Vectors[i].Target));
      SaveStream.Write(alph[i], sizeof(alph[i]));
      for j := 0 to n_Degree - 1 do
        SaveStream.Write(Vectors[i].Items[j], sizeof(Vectors[i].Items[j]));
    end;
    for i := 0 to n_Degree - 1 do
      SaveStream.Write(w[i], sizeof(w[i]));
    SaveStream.SaveToFile(FileName);
    Result := true;
  except
    Result := false;
  end;
  SaveStream.Free;
end;

function TSVM.LoadFromFile(FileName: string): boolean;
var
  i, j: longint;
  ReadStream: TMemoryStream;
begin
  ReadStream := TMemoryStream.Create;
  try
    ReadStream.LoadFromFile(FileName);
    ReadStream.Seek(0, soFromBeginning);
    ReadStream.Read(end_support_i, sizeof(end_support_i));
    N := end_support_i;
    setlength(Vectors, N); setlength(alph, N);
    ReadStream.Read(n_Degree, sizeof(n_Degree));
    for i := 0 to end_support_i - 1 do begin
      ReadStream.Read(Vectors[i].Target, sizeof(Vectors[i].Target));
      ReadStream.Read(alph[i], sizeof(alph[i]));
      setlength(Vectors[i].Items, n_Degree);
      for j := 0 to n_Degree - 1 do
        ReadStream.Read(Vectors[i].Items[j], sizeof(Vectors[i].Items[j]));
    end;
    for i := 0 to n_Degree - 1 do
      ReadStream.Read(w[i], sizeof(w[i]));
    Result := true;
  except
    Result := false;
  end;
  ReadStream.Free;
end;

destructor TSVM.Destroy;
begin
  setlength(Vectors, 0);
  setlength(alph, 0);
  setlength(w, 0);
  inherited;
end;

end.

---------------------------------------------------

调用窗体：

unit UnitMain;

interface

uses
  Windows, Messages, SysUtils, Variants, Classes, Graphics, Controls, Forms,
  Dialogs, StdCtrls, ExtCtrls, SMO_SVM;

type
  TFrmMain = class(TForm)
    ImgPoint: TImage;
    BtnTrain: TButton;
    BtnPredict: TButton;
    RdbRed: TRadioButton;
    RdbBlue: TRadioButton;
    procedure BtnPredictClick(Sender: TObject);
    procedure ImgPointMouseUp(Sender: TObject; Button: TMouseButton;
      Shift: TShiftState; X, Y: Integer);
    procedure FormCreate(Sender: TObject);
    procedure FormDestroy(Sender: TObject);
    procedure BtnTrainClick(Sender: TObject);
  private
    MySVM: TSVM;
  public
    { Public declarations }
  end;

var
  FrmMain: TFrmMain;

implementation

{$R *.dfm}

procedure TFrmMain.BtnPredictClick(Sender: TObject);
var
  i, j, X, Y: longint;
  point: TSingleExtendedArray;
begin
  setlength(Point, 2);
  for i := 0 to 320 do
    for j := 0 to 320 do begin
      point[0] := i / 320; point[1] := j / 320;
      if MySVM.Predict(Point) > 0 then
        ImgPoint.Canvas.Pixels[i, j] := clRed
      else ImgPoint.Canvas.Pixels[i, j] := clBlue;
    end;
  for i := 0 to MySVM.end_support_i - 1 do begin
    X := trunc(MySVM.Vectors[i].Items[0] * 320);
    Y := trunc(MySVM.Vectors[i].Items[1] * 320);
    if MySVM.Vectors[i].Target > 0 then begin
      ImgPoint.Canvas.Pen.Style := psClear;
      ImgPoint.Canvas.Brush.Color := $8080FF;
      ImgPoint.Canvas.Ellipse(X - 4, Y - 4, X + 4, Y + 4);
    end else begin
      ImgPoint.Canvas.Pen.Style := psClear;
      ImgPoint.Canvas.Brush.Color := $FF8080;
      ImgPoint.Canvas.Ellipse(X - 4, Y - 4, X + 4, Y + 4);
    end;
  end;
end;

procedure TFrmMain.ImgPointMouseUp(Sender: TObject; Button: TMouseButton;
  Shift: TShiftState; X, Y: Integer);
var
  point: TSingleExtendedArray;
  Outs: Extended;
begin
  setlength(Point, 2);
  point[0] := X / 320;
  point[1] := Y / 320;
  if rdbred.Checked then outs := 1 else outs := -1;
  MySVM.LearnExamples(point, outs);
  if rdbred.Checked then begin
    ImgPoint.Canvas.Pen.Style := psClear;
    ImgPoint.Canvas.Brush.Color := $8080FF;
    ImgPoint.Canvas.Ellipse(X - 3, Y - 3, X + 3, Y + 3);
  end else begin
    ImgPoint.Canvas.Pen.Style := psClear;
    ImgPoint.Canvas.Brush.Color := $FF8080;
    ImgPoint.Canvas.Ellipse(X - 3, Y - 3, X + 3, Y + 3);
  end;
end;

procedure TFrmMain.FormCreate(Sender: TObject);
begin
  MySVM := TSVM.Create;
  MySVM.Init(2, 1, 1, 2, 1, 8, klGsRBF);
  DoubleBuffered := True;
  ImgPoint.Canvas.FillRect(Rect(0, 0, 320, 320));
end;

procedure TFrmMain.FormDestroy(Sender: TObject);
begin
  MySVM.Free;
end;

procedure TFrmMain.BtnTrainClick(Sender: TObject);
begin
  MySVM.Learn;
  Application.MessageBox(PChar('There are/is ' + IntToStr(MySVM.Optimize) + ' support vector(s).'), 'Result');
end;

end.

-------------------------------

Related Url: http://research.microsoft.com/users/jplatt/smo.html